{'fraction_non_alphanumeric': 0.032942993179770945, 'fraction_numerical': 0.012675331360185305, 'mean_word_length': 4.876181474480151, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Jointly setting upper limits on multiple components of an anisotropic stochastic gravitational-wave background

With the increasing sensitivities of the gravitational wave (GW) detectors and more detectors joining the international network, the chances of detection of a stochastic GW background (SGWB) are progressively increasing. Different astrophysical and cosmological processes are likely to give rise to backgrounds with distinct spectral signatures and distributions on the sky. The observed SGWB will therefore be a superposition of these components. Hence, one of the first questions that will come up after the first detection of a SGWB will likely be about identifying the dominant components and their distributions on the sky. Both these questions were addressed separately in the literature, namely, how to separate components of isotropic backgrounds and how to probe the anisotropy of a single component. Here, we address the question of how to separate distinct anisotropic backgrounds with (sufficiently) different spectral shapes. We first obtain the combined Fisher information matrix from folded data using an efficient analysis pipeline PyStoch, which incorporates covariances between pixels and spectral indices. This is necessary for estimating the detection statistic and setting upper limits. However, based on a recent study, we ignore the pixel-to-pixel noise covariance that does not have a significant effect on the results at the present sensitivity levels of the detectors. We show that the joint analysis accurately separates and estimates backgrounds with different spectral shapes and different sky distributions with no major bias. This does come at the cost of increased variance. Thus making the joint upper limits safer, though less strict than the individual analysis. We finally set joint upper limits on the multicomponent anisotropic background using Advanced LIGO data taken up to the first half of the third observing run.


I. INTRODUCTION
A stochastic gravitational-wave background (SGWB) [1,2], created by the superposition of signals from unresolved astrophysical and cosmological sources [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17], are one among the most interesting targets for the current and upcoming gravitational wave detectors. Due to the nonuniform distribution of these astrophysical sources in the local Universe, the background can be significantly anisotropic [18][19][20][21][22][23][24]. Detection of such a SGWB seems promising in the near future and developing algorithms to probe the anisotropy in the background is, therefore, a highly relevant topic of research. Each source, generated by different phenomena, contributes to the anisotropic SGWB differently. Contributions from the individual components can be characterized by specific spectral shapes [25]. Upper limits have been placed on the anisotropic SGWB for primarily three spectral indices from the past observing runs of Advanced LIGO and Advanced Virgo [26][27][28][29].
These directional upper limits are placed for one spectral index at a time, neglecting the existence of the other components and thereby the covariance between differ-ent spectral components. The presence of one source contribution, however, can affect the estimation of the other components. Consider a situation where we filter the GW data for multiple SGWB components separately. We will end up overestimating the amplitudes of the individual components and underestimating the error bars. Hence, for the possibility of coexistence of multiple spectral components, one must estimate the amplitudes of individual components jointly, a joint-index multicomponent estimation is therefore necessary. This need does not supersede the existing results for single indices if one component is much stronger than the others. However, for scenarios where the strength of at least one component is not negligible, a joint estimation will be necessary. For instance, if we want to set an upper limit on the average eccentricity of extragalactic pulsars by observing the stochastic background created by them and, since we know that the same galaxy clusters also host the binary mergers (though the proportions may vary across the clusters), the latter cannot be ignored while setting the former upper limit [30]. On the other hand, one can also get limits on an unmodeled spectrum [29,31] by estimating the background in multiple independent frequency bins. This may not separate different signals but could reduce bias in the recovery that one might see if one uses a single spectrum over the whole analysis, though the variance in each bin will be significantly higher than that in the full broadband analysis.
Methods have been proposed in the past to disentangle arXiv:2106.09593v2 [gr-qc] 16 Nov 2021 these components [32][33][34][35]. These methods are appropriate for estimating the parameters associated with different SGWB models. However, with the increase in SGWB components, we have to deal with likelihood functions in multidimensional parameter space. In such situations, techniques like Markov-chain Monte Carlo or nested sampling will be needed to estimate the model parameters.
Recently, a new method has been proposed in Parida et al. [36], where the Maximum-Likelihood estimation problem for a multicomponent isotropic background was transformed to a linear deconvolution problem. This makes the estimation of the SGWB components trivial and computationally cheap. This novel component separation method was extended to the anisotropic SGWB case [37], where one can jointly estimate the amplitudes of the many components and set a directional upper limit. In this present work, we make use of the folded data [38] and the efficient PyStoch [31,39] pipeline to produce the first simultaneous 95% confidence level directional upper limits on multicomponent anisotropic SGWB from the Advanced LIGO's observational data taken up to the first half of the third observing run.
The paper is organized as follows: In Section II we discuss the component separation method. Sec. III briefly summarizes the dataset used. We describe the injection study performed to test the algorithm in Sec. IV. We also run the component separation algorithm on Advanced LIGO's observational data taken up to the first half of the third observing run. The details on this run along with the conclusions are described in Sec. V

II. METHODS
Assuming the SGWB to be unpolarized, Gaussian and stationary, the quadratic expectation value of the GW strain distribution h A (f,Ω) across different sky directions and frequencies can be written as where P(f,Ω) denotes the true anisotropy in the sky for a polarization A (from here on-wards, we assume that both polarizations contribute equally). Now we can write the SGWB energy density having the units of the dimensionless energy density parameter per steradian as It is also useful to write the above equation in terms of energy flux, which has units of erg cm −2 s −1 Hz −1 as Assuming that the source frequency power spectrum remains constant, and separable from its angular dependency, we can write the true sky as where H(f ) is the spectral shape of the background. If we assume that the SGWB is constituted by multiple components, denoted by the spectral index α, with a known spectral shape (we assume power-law spectrum throughout this work), here f ref is a reference frequency at which the detector is most sensitive is set as 25 Hz. Then the amplitude of the SGWB intensity can be rewritten as This expression differs from the usual form, P(f,Ω) = H(f ) P(Ω) used in GW radiometer algebra [38,40,41]. The corresponding changes in the algebra are given below.
We start with the usual cross-spectral density (CSD), C I (t; f ) := s * I1 (t; f ) s I2 (t; f ), of data from baseline I formed by the pair of detectors I 1 and I 2 , where s I1,2 are the short-term Fourier transforms (SFTs) of timeseries data obtained from the data segments of duration τ at time t. The usual expectation of the CSD [38] then gets modified to, where the sky-maps have been expanded in the basis e u (Ω) as P α (Ω) = u P α u e u (Ω) and γ I f t,u is the generalized overlap reduction function (ORF) [42,43], that accounts for the mismatch between the response functions of the detectors and the delay in signal arrival times, defined as, (8) In the above equation the polarizations are denoted by A = +, ×, and F A I1,2 (Ω, t) denotes the respective antenna pattern functions. Here ∆x I (t) is the separation vector between the detectors. In pixel basis, the ORF takes the form, whereΩ u is the direction of pixel u. The folded CSDs [38] average the SFTs over many sidereal days and can be assumed to be Gaussian. 1 Hence, for a signal model that allows for a superposition of multiple GWB components [as in Eq. (6)], and assuming that the detector noise spectra are well estimated, the likelihood function can be written as where the covariance of C I for different times and frequencies is given by the matrix N . This covariance matrix turns out to be diagonal, that is, data segments from different times and frequencies are uncorrelated. Data from different baselines I are also uncorrelated [44,45].
Then following the standard maximum-likelihood (ML) method for mapping the GW intensity [38,40,46] in a general basis, the clean map is given bŷ where, the "dirty map," and the Fisher information matrix, In the above equations P I1,2 (t; f ) are the one-sided noise power spectral density of the individual detectors. From Eq. (13), it is evident that the Fisher matrix is proportional to the noise-weighted inner product of the spectral shapes H α (f ) and H β (f ). The clean map estimator takes the simple form given in Eq. (11), when the spectral shapes are not degenerate. However direct inversion of the Fisher matrix could lead to numerical errors. This may demand some form of conditioning to perform the matrix inversion [46][47][48].

III. DATA
To perform component separation and to compute the upper limits using the joint-index formalism, we use the strain data 2 from the first (O1), second (O2), and firsthalf of third (O3a) observing runs of Advanced LIGO detectors located in Hanford (H) and Livingston (L) [49]. The fetched data have been processed and conditioned in the same way as was done in the previous analyses [26][27][28]45] to obtain the cross-spectral density and the corresponding variance, with a coarse-grained resolution of 1/32 Hz. Following the steps described in Ain et al. [38], these datasets are then folded into one sidereal day. We then perform the radiometer analysis "blindly" on the folded dataset to obtain the dirty map and the Fisher matrix corresponding to each observational run using PyStoch. One can obtain the combined Fisher matrix and dirty map from these data by following [2],

IV. INJECTION STUDY
We now demonstrate the component separation method for an anisotropic background using injections. In all the previous stochastic analyses performed on the LIGO's observing runs, three values of α have been considered. These spectral indices are related to different SGWB models: α = 0, corresponding to different cosmological models [50]; α = 2/3 is related to an astrophysical background which is dominated by compact binary coalescence [51]; and α = 3, indicating a flat strain spectrum [52] which stand for sources like millisecond pulsars and magnetars. Throughout this work, we will consider only these three spectral dependencies to estimate the SGWB amplitudes. Now for the three spectral index case, the clean map [Eq. (11)] can be obtained by solving the convolution equation in the matrix form The 3 × 3 matrix in the above equation denotes the coupling matrix C αβ uu formed by combining the Fisher information matrices for all possible combinations of the spectral indices. Here {α, β} ∈ {α1, α2, α3} denotes the spectral index which assumes different values for different contributing sources (each elements in the above equation represents either a vector or tensor with dimension of u and u ). We can now rewrite Eq. (11), the ML solution of the convolution equation, in terms of the coupling matrix asP However, in reality, the coupling matrix C is sparse in nature due to the covariance between the spectral shapes (these effects could further boosted by the pixel-to-pixel covarinace). Hence the direct inversion of C will lead to numerical errors. We resolve this by invoking a preconditioning of this coupling matrix. We achieve this by defining, where the C becomes a block matrix with a unit diagonal and the off-diagonal components are positive numbers.
To show the preconditioned coupling matrix elements, we chose one particular direction (100th pixel in a HEALPix resolution of N side = 16 following the RING [53] ordering scheme.), and then for that particular pixel C equates to As it is evident from the above coupling matrix example, the α = 0 and α = 2/3 components are highly correlated. Now we can rewrite Eq. (16), in terms of preconditioned coupling matrix aŝ where X β u = X β u / Γ ββ u u . It was recently shown [48] that at the current sensitivity level of the detectors, the contribution from the off-diagonal components of the pixel-to-pixel covariance matrix on the upper limits can be ignored. Following which we limit our analysis to retaining only the diagonal components of the pixel-to-pixel covariance matrix (however, as the detector sensitivity and number of detectors increase, we may have to incorporate these effects as well). Given we are neglecting the pixel-to-pixel covariance, Eq. (19) takes a simple form for a given pixel asP We can now solve the above linear convolution equations [Eq. (20)] with C to get the estimate of the SGWB sky.
To demonstrate our method, we first consider the O1+O2+O3a data from the HL baseline. Figure 1 shows the coupling matrix for all the three spectral indices (0, 2/3, 3) computed from this dataset. Now we consider source signals from different input maps to verify the joint-index multicomponent separation and hence the signal reconstruction. We use three types of toy models for injections: one containing a dipolelike sky pattern with α = 0; an extended source that mimics the diffuse part of the Milky Way Galaxy, whose spectral dependency is characterized by α = 2/3; and a sky-patch following α = 3 spectral shape. We work at a HEALPix resolution N side = 16 throughout the map-making process. Figure 2 shows the injected sky maps for all three cases along with the sky map obtained by combining these individual sources. The injection strengths we choose are such that the background signal is well above the noise level.
We now apply the usual radiometer map-making method [31] to the folded dataset with injections. In the radiometer analysis, we apply a time-varying phase delay to perform earth-rotation synthesis imaging and assign weights to the frequencies proportional to the expected source spectrum and inverse proportional to the detector noise spectrum [44], to obtain the dirty map. The dirty maps obtained for each choice of the spectral functions are shown in the first row of Fig. 3. After the dirty map creation, following the typical LIGO-Virgo-KAGRA (LVK) collaboration stochastic analyses, we neglect the pixel-to-pixel correlations 3 and obtain the semi-clean map from both the conventional ("single-index") and the new joint-index multicomponent ("joint-index") analysis.
The second row of Fig. 3 shows the semiclean SGWB sky maps for all the three spectral indices for the singleindex results. Here we assume that each SGWB spectral components are independent. Comparing with the injected source strength (see Fig. 2), it is apparent that we overestimate the background of each SGWB component for the single-index analysis. This is not surprising since we neglect the correlated background coming from other components when we choose a particular spectral function H(f ). Hence this severely biased method is not the correct way to extract the SGWB when the background components for the other spectral indices cannot be ignored.
We now use our proposed prescription and obtain the semiclean map by properly taking into account the covariance between different spectral shapes. For this, we solve the equation given in Eq. (20) by applying preconditioning [Eq. (17)] on the coupling matrix. For each pixel, we invert the 3 × 3 covariance matrix comparing different spectral indices. The third row of Fig. 3 shows thus obtained semiclean maps. It is evident from the maps that joint-index successfully undo the correlation between the three spectral shapes. To evaluate the accuracy of this recovery we have to use a certain quantitative measure, which can be challenging to define [40,47].
In this work we use the normalized mean squared error (NMSE) to quantify the deviation of the reconstructed semiclean map from the input map. If A, B are respectively the input map and the reconstructed map, NMSE is defined as, A better recovery is indicated by the value of the NMSE being closer to zero. Table I summarizes and compares the quality of the reconstructed sky map as a function of NMSE. The values of NMSE for the joint-index multicomponent separation case being close to 0, corroborates the importance of covariance between different components.
To prove that the NMSE behavior we showed in Table I is not obtaining by chance, we repeat the procedures described in the injection study for 1000 noise realizations. We then compute the NMSE for each case. The histograms obtained from this study employing both single-index and joint-index estimations for all the three spectral indices are shown in Figure 4. It is evident from the histograms that, for all these cases, NMSE for jointindex estimation is close to zero indicating a better source reconstruction. On the other hand, the single-index estimations are biased, indicated by the scattered nature of NMSE which is always significantly greater than the one from joint-index estimation. This further confirms the better source recovery of the joint-index multicomponent estimations.
Given the current sensitivity levels of the detectors, when we analyze real data, we will be more interested in setting an upper limit on the SGWB components. Here, we investigate the effect of the joint-index multicomponent separation method on the upper limit. We first consider three spectral shapes and three sky patterns as in Fig. 2 and use the same injection strength. We then compute the 95% confidence upper limit [45,54] on the amplitude of the recovered signals from both methods. Fig. 5 show the results obtained from this study. In the top panel, the estimated amplitudes from both singleindex and joint-index estimations are shown. Comparing with the injection strength, single-index estimations always overestimates the amplitude of the individual components. We explore further to see the implication of a biased amplitude estimation on the upper limit. We first compute the variance of single and joint-index estimators. As we noted earlier, the variance or the error bar is underestimated for the single-index analysis. From earlier studies [36,46,48], it is shown that when the parameters of the analysis increase, the variance of the corresponding estimator also increases (see Appendix A for more details). Here, for the joint-index analysis, we observe the same behavior. These responses are depicted with the aid of histograms in the middle panel of Fig. 5. We then compute the 95% confidence upper limit at every pixel on the sky using both methods. The results are shown in the bottom panel of Fig. 5. The histogram in the bottom panel shows the difference between the upper limit we computed and the injection strength. However, if we consider, α = 2/3 case in our injection study, we can see that a considerable amount (∼ 11%) of the points in the histogram fall under the negative region. This indicates that, even though the single-index estimate appears to provide a stronger upper limit, it violates the requirement for a 95% confidence upper limit, at least in certain generic cases, where for more than 5% pixels (∼ 11% here) the upper limit is lower than the injected value (To be more rigorous, one should create several signals and noise realizations and find out in what fraction of cases the upper limits are lower than the injected value. This will however be computationally expensive. If one assumes that most pixels are reasonably independent, one would expect 95% of the pixels to have upper limits higher than the injected value). This exercise suggests that even if the detectors are not sensitive enough to detect SGWB, the joint-index multicomponent estimator provides safer upper limits when one cannot ignore the existence of more than one component.

V. RESULTS AND CONCLUSION
We now apply the joint-index multicomponent separation method to the folded dataset derived using the observational data from aLIGO's H and L detectors, taken up to the first half of the third observing run. Similar to what we did for the injection study, we first compute the dirty maps from the O1+O2+O3a observation data by employing the optimal filtering [44]. We then construct the coupling matrix by assuming three spectral indices α = 0, 2/3, 3 from the same dataset. To obtain the true SGWB sky, we solve the matrix equation given in Eq. (15). Here we account for spectral index covariance.
Since the noise level from the observed data is still larger for signal detection, the 95% upper limits are placed on the gravitational wave energy flux (units of erg cm −2 s −1 Hz −1 ) [26][27][28]45]. To further demonstrate the importance of a joint-index estimation, we choose four pixels corresponding to four random sky directions and explore the effect of the covariance between the spectral indices. We compute both single and multicomponent upper limits for these directions (given the dataset is dominated by noise and the existence of pixel-to-pixel covariance, it is not easy to show the corner plot for every pixel). The one-dimensional marginalized probability distribution and the error ellipses are shown in the corner plot Fig. 6. These corner plots show that the upper limit gets underestimated for different pixels in the sky if we rely on single-index analysis. We can also see that this deviation, compared to the joint-index upper limit, is substantial in some cases (see top or bottom right corner plot). From this figure, it is evident that the effect of spectral covariance on the estimation of SGWB components is significant. Given the results presented here, in order to account for these effects, it demands the use of joint-index multicomponent estimation in the SGWB search.
In Fig. 7, we provide the 95% upper limit range for the gravitational wave energy flux for all the sky directions. The top panel shows the upper limit sky map obtained using the single-index estimation. Whereas the joint-index multicomponent estimation is shown in the bottom panel. All these sky maps are consistent 4 with the upper limit range quoted in the most recent directional limits from the LVK collaboration [28]. We also showed the 95% credible upper limit sky maps by simultaneously considering two spectral indices in Figure 8. Here each column represents the results from all three combinations of spectral indices. We summarize these results in Table II for easy comparison. Here the upper limits on individual source components are placed following the single-index estimation, the SGWB amplitudes FIG. 6. The estimates of three SGWB components (for α = 0, 2/3, 3) using O1+O2+O3a folded dataset for four random sky directions (pixel number = [71,200,1566,1999] in a HEALPix resolution of N side = 16 following the RING [53] ordering scheme) are shown. The purple (dashed) 2σ error bar shows the result from a single-index analysis, whereas the cyan (solid) one represents the error bar from the joint-index analysis. The shaded and unshaded regions represent the 1σ and 2σ confidence intervals respectively. The one-dimensional plot depicts the marginalized probability distribution for each component. Here the solid vertical line denotes the mean, and the dotted vertical line represents the 95% confidence upper limit after marginalizing over the calibration uncertainty [54] and assuming a uniform prior on SGWB flux. As evident from Figure 5, when there is an injection, the intersection of cyan (solid) error bars pass through the point marking the injection value. The plot establishes that the effect of covariance between the spectral shape on the estimates are significant and hence the need for a multicomponent analysis.
are overestimated and hence biased. Therefore, the unbiased results from the joint-index estimation give a more accurate representation of the SGWB background upper limit from individual source categories if non-negligible multiple components are present together.
Relating the SGWB signal to the sources that contribute to it is a challenging problem. While extracting the individual source contribution, minimizing the bias in the estimators is important. With the increase in detector sensitivity, an unbiased estimator, like the one shown in this paper will act as an essential tool. The method and results shown in this work will set a benchmark for this future problem, which is likely to become an essential task after the first detection of SGWB.
This analysis can be extended to multiple directions. Since the now-standard analysis pipeline PyStoch can  Table II FIG. 8. Sky map representing 95% upper limit on the SGWB energy flux (erg cm −2 s −1 Hz −1 ) obtained from O1+O2+O3a folded dataset by jointly considering two spectral index. From left to right each column represent upper limit sky maps from joint-index analysis considering (0,2/3), (0,3), and (2/3,3) spectral index combinations. All the maps are represented as a color bar plot on a Mollweide projection of the sky in ecliptic coordinates. For the visualization purpose we smooth the upper limit sky maps with a 10 degree Gaussian beam and adjusted the color bar scale using the values from Table II. perform a model-independent multi-baseline search for the anisotropic SGWB straightforwardly, it will be interesting to combine it with the proposed method. This will enable us to simultaneously estimate the spectral, frequency, and angular dependence of any SGWB signal. This will be the keystone to future component separation efforts. It will be also interesting to look at this component separation problem using the spherical harmonic (SpH) decomposition techniques. Given the information about the coupling of different SpH modes are encoded in the SpH fisher information matrix, it could potentially help us in separating the different spectral shape contributions to the SGWB backgrounds from different multipoles. At this stage, since the LVK stochastic analysis set the upper limit assuming three source categories, characterized by three spectral shapes, the component separation method presented in this paper will augment the existing analysis.  Figure 7 and 8. The increase in upper limit exhibited by the joint-index estimation is due to the increase in the variance of the corresponding estimator.

ACKNOWLEDGMENTS
The authors thank Patrick Meyers for carefully reading the manuscript and providing valuable comments. This work significantly benefited from the interactions with the Stochastic Working Group of the LIGO-Virgo-KAGRA Scientific Collaboration. We acknowledge the use of IUCAA and Caltech LDAS clusters for the computational/numerical work. J. S. acknowledges the support by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant No. JP17H06361. S. M. acknowledges support from the Department of Science and Technology (DST), India, provided under the Swarna Jayanti Fellowships scheme. This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation. This research has also made use of data obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org/ ), a service of LIGO Laboratory, the LIGO Scientific Collaboration, and the Virgo Collaboration. LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN), and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. This paper is assigned the LIGO document control number LIGO-P2100169. Parts of the results in this work make use of the colormaps from the CMasher package [55] Appendix A: Number of components and variance Increasing the number of components increases the uncertainty or variance associated with each component estimations. It has been shown that [36] if we neglect the covariance between spectral shapes, we underestimate the error bar. On the other hand, the error bar increases when we simultaneously estimate the individual spectral shapes. We have shown this behavior in Figure 9, where we used the joint-index analysis by considering only twocomponent simultaneously. By comparing the histogram of the variance with the one for three-component case (see Figure 5) we can see that the variance reduces for all the two-component cases (it is also worth noting that, here also, the single-index analysis for α = 2/3, many points in the upper limit histogram (second panel from top) fall under the negative region). This suggests that there is a trade-off between accurately extracting the individual components and minimizing the uncertainties. In general, the addition of detectors to the network improves this trade-off problem. However, given the current sensitivity of Virgo detector [56], adding it to our analysis will not contribute considerably. To demonstrate these effects, we consider all the three baselines (HL, HV, LV) and eigenvalues of their respective coupling matrix (here we consider both pixel-to-pixel and spectral index covariance). Eigenvalues obtained from the singular value decomposition are plotted and compared with the combined baseline result (HLV) in Figure 10. It is evident from this figure that, given the sensitivity of the Virgo detector, the most dominant contributions to the eigenvalues are coming from the HL baseline. This situation can improve drastically when detectors like KAGRA [57,58] and LIGO-India [59,60] start their observation run.