Stochastic gravitational wave background mapmaking using regularised deconvolution
SStochastic gravitational wave background mapmaking using regularised deconvolution
Sambit Panda, ∗ Swetha Bhagwat,
2, 3, † Jishnu Suresh,
4, 5, ‡ and Sanjit Mitra § Birla Institue of Technology and Science (BITS) Pilani, Pilani, Rajasthan 333031, India Dipartimento di Fisica, Sapienza Universita di Roma, Piazzale Aldo Moro 5, 00185, Roma, Italy Syracuse University, Syracuse, NY 13244, USA Institute for Cosmic Ray Research (ICRR), KAGRA Observatory,The University of Tokyo, Kashiwa City, Chiba 277-8582, Japan Inter-University Centre for Astronomy and Astrophysics (IUCAA), Pune 411007, India
Obtaining a faithful source intensity distribution map of the sky from noisy data demands incorpo-rating known information of the expected signal, especially when the signal is weak compared to thenoise. We introduce a widely used procedure to incorporate these priors through a Bayesian regular-isation scheme in the context of map-making of the anisotropic stochastic GW background (SGWB).Specifically, we implement the quadratic form of regularizing function with varying strength of regu-larization and study its effect on image restoration for different types of the injected source intensitydistribution in simulated LIGO data. We find that regularization significantly enhances the qualityof reconstruction, especially when the intensity of the source is weak, and dramatically improvesthe stability of deconvolution. We further study the quality of reconstruction as a function of reg-ularization constant. While in principle this constant is dependent on the data set, we show thatthe deconvolution process is robust against the choice of the constant, as long as it is chosen froma broad range of values obtained by the method presented here.
PACS numbers: 04.80.Nn, 95.55.Ym, 98.70.Vc
I. INTRODUCTION
A stochastic gravitational wave background(SGWB) [1, 2] can arise from the gravitationalwaves (GW) produced by unresolved astrophysical andcosmological sources. For instance, a large number ofunresolved distant compact binary coalescences [3]and millisecond pulsars in galaxy clusters [4–7] canproduce an SGWB. The background can be significantlyanisotropic due to the non-uniform distribution ofastrophysical sources in the local universe [8]. Basedon the rates estimated from the observed GW signals,the detection of SGWB created by compact binariesseems promising in the near future [3]. With thedetectors becoming progressively sensitive, along withnew detectors coming online and with the promise ofmultiple next-generation detectors, making the map ofthe SGWB-sky may soon become a reality. The detec-tion of an anisotropic SGWB will open an independentwindow to the universe and enable us to probe persistentGW sources which are electromagnetically dark.Several algorithms have been proposed in the past toprobe both isotropic and anisotropic SGWB [9–14]. Herewe focus on the directed radiometer search [15–17]—astandard method to construct a sky-map of SGWB in apixel-basis from the output of a network of ground-basedinterferometric detectors. A raw sky-map produced bythis analysis is referred to as the ‘dirty map’. A dirty ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] map is a convolution of the true source map with theantenna response function of the detectors and containsadditive noise. To obtain an estimate of the true sourcedistribution in the sky, known as ‘clean map’ one needsto undo the effects introduced by the detector’s antennapattern in the dirty map through a procedure called de-convolution or ‘cleaning’. However, the process of de-convolution, although mathematically well defined, neednot be numerically stable and the procedure itself mayintroduce an additional source of error. In the case ofSGWB mapmaking, it is known that the deconvolutionprocess introduces a considerable amount of noise in thecleaned map because of insensitive modes of the beam.Consequently, it has been a common practice to use thedirty map directly to obtain scientific results. Here wepropose a scheme to efficiently obtain clean maps using‘regularised’ deconvolution by incorporating prior knowl-edge of the signal in a Bayesian framework [18–20].Since the source map is unknown a priori , one needsa model independent prescription for implementing theregularisation scheme into the mapmaking procedure.Our method uses broad features of the expected char-acteristics of the source and the noise as priors. Incorpo-rating the priors in this specific context is algebraicallyequivalent to adding a ‘regularisation function’ to the log-Likelihood. Motivated by earlier successful attempts [18–20], we use quadratic forms for the regularisation func-tion. In particular, the regularization function we imple-ment are constructed such that they penalize the ‘norm’or the ‘gradient’ of the cleaned map. We demonstratemapmaking using this procedure for two specific cases -a point source and an extended source closely mimickingthe diffuse part of the Milky Way galaxy. In addition,one needs to make a choice on how much the clean map a r X i v : . [ g r- q c ] J un is allowed to depend on prior knowledge. This choice iscoded in a quantity called ‘regularization constant’, thatbalances the credence on the data and the prior. We alsodemonstrate that our procedure is not very sensitive tothe exact value of regularization constant, provided it ispicked from a range decided by following the prescriptionpresented in this paper.The paper is organized as follows. Section II briefly re-views the GW radiometer analysis. Section III providesa detailed description of sky map reconstruction usingregularized deconvolution. In section IV, we discuss thenumerical implementation for in our study. In Section V,we present the results and show that the regularisationprocedure improves the quality of deconvolution statis-tically. In section VI, we discuss the implications of ourmethod on SGWB mapmaking and its immediate appli-cability to the data from the current detectors. II. GW RADIOMETER ANALYSIS
The GW radiometer analysis [15–17, 21] is currentlythe standard technique to probe an anisotropic SGWBproduced by unresolved astrophysical and cosmologicalsources. Its fundamental principle is based on the as-sumption that noise in geographically well-separated de- tectors are uncorrelated. Therefore, one can use thecross-correlation between detector outputs as a statisticalmeasure of the stochastic GW signals. We use the cross-correlation statistic obtained by applying a direction de-pendent filter that accounts for the phase delay in thearrival of the GW signal at the detector locations, whichis adjusted as the earth rotates. The statistic is averagedover the time segments with inverse noise weights to im-prove the signal-to-noise ratio (SNR) [16, 17]. The filteris optimized on small time segments, typically 32 −
192 seclong, and it also depends on the expected power spectraldensity (PSD) of the source and the detector noise. Fur-ther, it is generally assumed that each point source hasthe same PSD, and strengths in the two GW polariza-tions are equal. Then, the resultant dirty map can bethought of as a convolution of the true sky map with aneffective kernel, often called as the beam function ( B ). Itcan be expressed as, D ( ˆΩ) = (cid:90) d ( ˆΩ (cid:48) ) B ( ˆΩ , ˆΩ (cid:48) ) S ( ˆΩ (cid:48) ) + n ( ˆΩ) , (1)where S is the source map and n is the noise. The beamcaptures the effect in the direction ˆΩ due to a source inthe direction ˆΩ (cid:48) and can be written as, B ( ˆΩ , ˆΩ (cid:48) ) = Λ( ˆΩ) (cid:90) dt (cid:90) −∞∞ df H ( f ) P ( t ; f ) P ( t ; f ) (cid:88) A (cid:48) F A (cid:48) ( ˆΩ (cid:48) , t ) F A (cid:48) ( ˆΩ (cid:48) , t ) (cid:88) A F A ( ˆΩ , t ) F A ( ˆΩ , t ) e πif ˆΩ · ∆ x ( t ) c , (2)where, Λ is a normalization constant, H ( f ) is the sourcePSD, P , ( t ; f ) are the one-sided PSD of detector noiseand ∆ x ( t ) = x ( t ) − x ( t ) is the geometrical separationbetween the detectors located at x ( t ) and x ( t ) respec-tively. The polarizations are denoted by A = + , × and F A , ( ˆΩ , t ) are the antenna functions of the detectors.In practice, the convolution is done in a discrete basis.An extended source is modeled as a smooth collectionof point sources, i.e., a collection of bright pixels. Themaps are thus represented by vectors living in an n pix -dimensional space, where n pix is the number of pixel inthe map, and the beam becomes a matrix. Further, onecan normalise the map in such a way [17] that the beammatrix becomes symmetric and the pixel-to-pixel noisecovariance matrix becomes proportional to the beam ma-trix, which makes it convenient for numerical manipula-tions. The pixelised dirty map D can then be cast as aset of linear convolution equations, D = B · S + n , (3)where S is the pixelised source map and n represents(convolved) detector noise. The i th row of the matrix B is the response function of the radiometer when pointed to the i th pixel in the sky. Obtaining an optimal estimateof S , the ‘clean’ map ( ˆ S ), from the above convolutionequation is the primary goal of this paper.A graphical representation of a typical beam matrixfor the baseline created by the two LIGO detectors withthe sky tessellated in n pix = 3072 pixels is depicted inFig. 1; the colour bar indicating the values of the matrixelements. While the matrix is diagonally dominated, in-dicating that the maximum contribution in a map comesfrom the corresponding pixel in the source map, there isa significant number of off-diagonal elements which arisebecause the beam is broad.The pixel-to-pixel noise covariance matrix of the dirtymap has a non-trivial form. With the choice of normal-isation mentioned above, the noise covariance matrix N becomes proportional to B and can be written as, N = κ B , (4)where κ is a constant that depends on the exact nor-malisation convention being used, which can involve the FIG. 1: A characteristic beam matrix for the LIGO Hanford-Livingston baseline at ∼ ◦ pixel size (HEALPix a n nside =16). Each row of the matrix is the antenna response functionfor the pointing direction. Since the GW radiometer receivesmaximum contribution from the pointing direction, this ma-trix is dominated by the diagonal elements. Here the stripesare related to the isoLatitude pixelization scheme - the in-dices of the neighboring pixels at different latitudes differ bythe total number of pixels on that latitude. The beam alsotakes negative values. a http://healpix.sf.net duration of each time segment and frequency bin size. III. REGULARIZED DECONVOLUTION
To obtain an estimate for the true skymap, ˆ S , the di-rect inversion of the convolution equation,ˆ S = B − D , (5)boosts the noise associated with the low sensitivity modesof the beam, resulting in a sky reconstruction dominatedby numerical noise. This can be understood by lookingat the very low singular values of the matrix B by per-forming a singular value decomposition (SVD) [17]. Thisleads to a small condition number (that is, the beam isill-conditioned), rendering the direct inversion unreliable.Therefore, given a D , there may not be a unique numer-ical solution ˆ S . The weak intensity of the SGWB signalwith respect to the detector noise level further enhancesthis issue. As a result, a direct inversion of the convolu-tion equation does not lead to a reliable solution. The qualitative features of our analysis is independent of theexact value of the proportionality constant κ . The optimum reg-ularisation constant, which is numerically determined, absorbsthis factor. However, problems of this nature exist in many fieldsand a traditional way to tackle this has been through in-corporating prior information about the system. For in-stance, a similar discussion on image reconstruction couldbe found in the context of strong gravitational lensing inSuyu et al. [18]. In a Bayesian framework, one can imple-ment various regularization schemes as priors and incor-porating them in the deconvolution routine is expectedto remedy the ill-posed nature of Eq. (5).
A. Most Likely Solution : S ml In a Bayesian framework, the problem can be formu-lated as finding the most likely solution. Assuming aGaussian noise model the likelihood function is given by, P ( D|S , B ) = 1 Z L e − ( D−BS ) T N − ( D−BS ) . (6)Here Z L is the normalization constant. The signal thatmaximizes the above likelihood function gives an esti-mate ˆ S , which is called the most likely solution S ml .Another perspective to this problem can be obtainedby looking at the χ . Using the standard χ estimator,the likelihood function can be written as, P ( D|S , B ) = 1 Z L e − χ , (7)and S ml can be interpreted as the solution that minimizesthe χ , ∇ χ ( S ml ) = 0 . (8)Solving the above equation one gets, S ml = ( B T N − B ) − ( B T N − D ) . (9)Substitution of N = κ B then reduces the above formulato Eq. (5). Thus, likelihood maximization in the contextof SGWB map-making reduces to an ill-posed problem.This is because the set of linear equations that describedeconvolution are not independent, i.e., the number offree parameters in the problem is much more than thenumber of constraint equations. This invokes the risk ofover-fitting by minimizing χ to an unrealistically smallvalue. One has to incorporate regularization for a reli-able minimization of χ . Regularization serves as a setof additional constraints to the set of linear equations,Eq. (3), and prevents over-fitting of data. B. Most Probable Solution : S mp In a Bayesian framework, the prior distribution func-tion incorporates additional information on the signal inthe form of a function R ( S ) which we refer to as theregularization function. Given R ( S ) and a regularizationstrength µ , the prior distribution can be written as, P ( S| R ( S ) , µ ) = 1 Z P e − µR ( S ) , (10)where Z P is the normalization for the probability. UsingBayes theorem, the posterior probability of getting thesignal S , given the dirty map D , the beam B and theregularization function R ( S ) is as follows, P ( S|D , B , µ, R ( S )) = P ( D|S , B ) P ( S| R ( S ) , µ ) P ( D| µ, B , R ( S )) . (11)Unlike P ( D|S , B ), which depends only on data D , P ( S|D , B , µ, R ( S )) evaluates the probability of solution S by combining the information in both D and R ( S ),through the posterior distribution, P ( S|D , B , µ, R ( S )) = 1 Z P e − M ( S ) . (12)where Z P is the Bayesian evidence. M ( S ) is defined as, M ( S ) := χ ( S ) + µR ( S ) , (13)which contains two competing function χ and R ( S ). S mp is obtained by maximizing the posterior probability,that minimizes M ( S ), the linear combination of thesecompeting functions, via solving, ∇ M ( S mp ) = ∇ ( χ ( S mp ) + µR ( S mp )) = 0 . (14)This is analogous to Eq. (9). Let λ := µ ( κ ∆ t/
4) and C be the Hessian of the regularization function R ( S ), thenit can be shown that S mp acquires the following simpleform [c. f. Appendix A for the detailed derivation], S mp = ( B + λ C ) − D . (15)From the above equations, regularization can also beviewed as modifying the convolution kernel such that theproblem becomes less ill-posed. In short, the reconstruc-tion of the source map without a prior corresponds toevaluating the expression for S ml , while regularized de-convolution corresponds to finding S mp . C. Regularization function
We choose to investigate the effect of norm and gra-dient regularization functions on reconstruction of thesource map.
1. Norm Regularization
Norm regularization introduces a preferential bias to-wards solutions that minimizes the norm of the map. This is seen to have a noise suppression effect, especiallyin the case of point-like sources. R norm ( S ) = 12 n pix (cid:88) i =1 S i . (16)It is clear that the Hessian matrix for this regularizationform is the identity matrix, C ij := ∇ S i ∇ S j R norm ( S ) = δ ij . (17)
2. Gradient Regularization
Gradient regularization incorporates a preference to-wards smooth source reconstruction by penalizing the in-tensity difference between the neighboring pixels. Thisis motivated by the idea that the spatial fluctuations inthe intensity of the noise will be more than the spatialfluctuation in the signal, especially if the signal has anextended pattern in the sky, that is, it prefers minimumvariation of intensities in the reconstructed map. Thegradient is minimized by using the following form of reg-ularization, R grad ( S ) = 12 n pix (cid:88) i =1 n i (cid:88) k =1 ( S i − S j ik ) , (18)where j ik is the pixel number of the k th nearest neighbourto the i th pixel and n i is the number of pixels touchingthe i th pixel, n i = 8 if each pixel has four sides. TheHessian for this form becomes, C ij := ∇ S i ∇ S j R grad ( S ) = δ ij − n i δ ij . (19)Notice that these two regularization functions do not con-tain any information about a specific source intensity dis-tribution and it is in this sense that they are uninforma-tive priors. These regularization schemes use knowledgeof statistical properties of both the noise and the ex-pected source distribution. However, we emphasize thatregularization introduces a non-zero bias. Derivation ofthe bias estimate is presented in Appendix B. D. The Strength of Regularization
The strength of regularization λ decides the relativeweights between the goodness-of-fit and the bias towardsour prior knowledge. It is crucial to choose a λ thatstrikes an optimal balance between the two. While verylow values of λ increase the risk of over-fitting the data,setting λ to large numbers will lead to a highly biasedsolution. In principle, one could systematically computethe optimal value of λ given a data set by maximizinglog P ( λ |D , B , R ( S )) with respect to log λ .In this section, we present a systematic prescription topick an optimal value of λ , which in principle depends onthe data itself.However, generally this happends to be a transcenden-tal equation that needs to be solved iteratively. Fur-thermore, it involves calculating S mp and is, therefore,computationally very expensive. Instead, for our study,we perform the deconvolution using a range of λ andpick the λ that maximizes the estimator of the qualityof reconstruction. [c.f the variation of the estimator asa function of λ as presented in Fig. 2 and 4]. However,we confirm that the λ that maximizes the estimator ofquality of reconstruction also satisfies Eq. ( ?? ) up to nu-merical tolerance. Furthermore, in Section V we arguethat our deconvolution procedure is robust to the choiceof λ as long as it is within a range which can be deter-mined by following the above procedure. E. Measure of quality of recovery
There is no unique prescription to construct a quan-titative estimator that measures the quality of a recon-structed map. If one considers the source map and theclean map as two vectors living in a vector space of di-mension n pix , a measure that quantifies the quality of thereconstructed map would be to calculate the separationbetween the two vectors. However, while measuring thisdistance, there is no unique prescription for choosing themetric on the vector space which contains these vectors.If we want to use distance as an estimator for detection ofa pattern in the sky, then a good choice of metric wouldbe one that minimizes false dismissals and false alarms.Construction of this estimator is beyond the scope of thisstudy and will be investigated in a future work.For our study, we construct an estimator called theNormalised Scalar Product, NSP, which quantifies the de-viation of the source map and the recovered map throughan inverse norm weighted Euclidean inner product. Wedefine NSP as, NSP = A · B (cid:112) (cid:107) A (cid:107)(cid:107) B (cid:107) , (20)where, A , B in our case are the source map and the recon-structed map respectively. We show that this quantifiesthe goodness of recovery sufficiently well, in the contextof this study for an extended source skymap. However,we caution the reader that it may be possible to constructbetter measures by a more careful choice of the metric indefining the inner product.Furthermore, we also use the Normalised MeanSquared Error (NMSE) as another independent measureto quantify the deviation of the reconstructed map fromthe source map. NMSE is mathematically defined as,NMSE = (cid:107) A − B (cid:107) (cid:107) A (cid:107) . (21)We find that NSP is a better quantifier for an extendedsource and NMSE is more suited for a point source. This is because NSP gives zero weight to pixels where no in-jections were made; hence noise in those pixels in therecovered maps are ignored, which may not be appropri-ate when a localised source is under consideration. Onthe other hand, when we are interested in an extendedpattern on the sky, matching of the recovered pattern ismore important than the overall normalisation, which iscaptured well by NSP. Needless to mention, better recov-ery is indicated by higher NSP or lower NMSE. IV. NUMERICAL IMPLEMENTATION
In our study, calculation of B is the most computa-tionally expensive step; it is equivalent of making onedirty map for each pixel by placing a unit point sourceat that pixel. For this we use the software pipelinecalled PyStoch [22] that utilizes advantages of foldeddata [23, 24] and the HEALPix [25] pixelization scheme.By considering the optimal resolution required for the ra-diometer analysis for the two LIGO detectors at Hanfordand Livingston, we choose n side = 16 in the HEALPixscheme. This corresponds to a pixel width of ∼ ◦ andthe whole sky is tessellated into n pix = 3072 pixels.Each injection is done such that every pixel of the mapis assigned a coefficient which corresponds to the inten-sity of that pixel in an arbitrary unit. Although arbitrary,this unit is consistent across all the plots presented in thispaper. Therefore, it can be used to compare the relativestrengths of the injected signals. Note that defining anSNR for these maps turns out to be tricky due to thenon-trivial pixel-to-pixel covariance of the noise.The time duration of each segment, ∆ T , is taken as52 sec. The upper cut-off frequency is set to 512 Hz.The source is assumed to have a flat PSD, H ( f ) = 1. Wefollow recipes given in Ain et al. [22] and Mitra et al. [17]to generate dirty maps with noise from injected sources.We implement Bayesian regularized deconvolution toreconstruct the source map from the dirty map. We varythe parameter λ over a range of 1 to 10 in logarithmicintervals and pick a value that (nearly) maximizes NSPor minimizes NMSE, as shown in Fig. 2 and 4 for differentinjection strengths. We later demonstrate that the qual-itative results are weakly sensitive to the choice of λ . Toperform the inversion of Eqs. (9) and (15) we use an in-built conjugate gradient solver ( c gs) in the Python SciPy package. We set a tolerance of 10 − for our study and amaximum iteration of 50 − V. RESULTS
We demonstrate the capabilities of our method for theextended sources as well as for the localized point-likesources with varying intensities. We apply both the gra-dient and the norm regularized deconvolution schemesfor each of these cases to obtain the clean maps. Here we
Source iter
No-reg λ NSP dirty
NSP
No-reg
NSP
Norm
NSP
Grad
NMSE dirty
NMSE
No-reg
NMSE
Norm
NMSE
Grad
Strong extended 12 5 0.9926 0.8437 0.8360 0.8798 0.0155 0.2887 0.3013 0.2259Weak extended 5 50 0.9033 0.6587 0.6448 0.8407 0.2486 0.6169 0.5851 0.2932Very weak extended 3 500 0.6678 0.3817 0.3533 0.7906 1.5539 1.3300 0.8752 0.3756Strong point 3 100 0.9384 0.2530 0.2544 0.2261 0.1447 0.9475 0.9409 0.9491Weak point 3 1000 0.5186 0.0711 0.0965 0.0820 3.6186 2.9303 1.0379 1.0249Very weak point 3 1000 0.3243 0.0381 0.0536 0.0547 14.4745 9.4343 1.2914 1.1830TABLE I: Quantitative measures of goodness of reconstruction in terms of NSP and NMSE of sky-maps shown in Figures 3 and 5.NSP (better measure for extended sources) and NMSE (better measure for point sources) are quoted for recovered maps obtainby no deconvolution (comparing dirty map to beam convolved injected map), unregularised deconvolution, and norm & gradientregularised deconvolution. The number of iterations for unregularised deconvolution (iter
No-reg ) and regularisation strength ( λ )are also listed. Except for strong sources, incorporating regularization significantly improves the quality of reconstruction. emphasize that gradient regularization is better tailoredfor an extended source distribution and the norm regu-larization scheme performs better for a point-like source.Furthermore, NSP is a better quantifier for the qualityof deconvolution for an extended source, while NMSEis better suited for a point-like source. However, irre-spective of the true source distribution, we present thevalues of both NSP and NMSE for quantifying the per-formance of the gradient and norm regularisation in Ta-ble I. For brevity, we include only selected plots in thepaper that correspond to (1) Gradient regularisation forextended sources with the NSP as the quality measureand (2) Norm regularisation for point sources with theNMSE as the quality measure. We then perform a simu-lation to show how the choice of λ without any fine-tuninggenerically improves the quality of deconvolution. A. Extended source: Gradient regularisation
Here we present the results of implementing gradientregularized deconvolution for extended sources. Since thepixel-to-pixel variation of noise is expected to be muchhigher than that of the source, one is motivated to choosea gradient regularization scheme for recovering extendedsource patterns (a preference to smoother reconstructedmap). The results are qualitatively presented in Fig. 3and quantitatively in Table I.Each column of Fig. 3 corresponds to an increasinglydiminishing injected signal strength, which can be seenin the magnitude in the colour bar of the top row of thefigure indicating the pixel intensities in an arbitrary unit.The second row shows the corresponding dirty maps,and the third row shows the unregularised clean maps.Unregularised deconvolution here tends to diverge if thenumber of iterations is increased beyond a certain value(the value depends on the injection). We plot the NSP ofthe recovered map against the number of iterations andchoose the number that corresponds to the best quality,i.e., maximum NSP. An example plot for a strong injec-tion is provided in the top left panel of Fig. 7. When thesignal is very weak, unregularised deconvolution starts diverging from the first iteration, making the reconstruc-tion unreliable. Nevertheless, we set the number of it-erations to 3 in such cases for plotting and comparisonpurposes.The fourth row of Fig. 3 shows the results of regulariseddeconvolution. We determine the optimal regularisationconstant λ from Fig. 2 based on the process described inSection IV. Weaker the source, stronger the regulariza-tion necessary. This is consistent with what one expectsnaively. We set λ = 5 , ,
500 respectively for the in-jection described here. For weak sources, regularisationintroduces significant improvement. The dirty and un-regularised clean maps for weak injections are dominatedby noise, while regularisation brings out some features ofthe injected source. Table I shows this quantitatively.For instance, we see that for a very weak injection NSPincreases from 0 .
38 to 0 .
79, while for a strong source,regularisation performs sub-optimally (sometimes wors-ens the result) and introduces a bias.
B. Point source: Norm regularisation
Norm regularization is motivated by the idea that alarge part of the sky generally does not contain anysource. The gradient regularization scheme is not op-timal for such localized sources, as it smears out sharpfeatures in the reconstructed maps. Here we present theresults obtained from norm regularized deconvolution ap-plied to point source injections of different intensities.The maps are presented in Fig. 5 and the quantitativeresults are provided in Table I.The first row of Fig. 5 shows the injected point sourcesfor three cases with diminishing strengths (from left toright). The second row shows the corresponding dirtymaps. The third row shows clean maps obtained by un-regularised deconvolution. Here the implementation ofconjugate gradient is unstable and the results get sat-urated by numerical noise for even a small number ofiterations. We therefore set the number of iteration to 3for all the cases here.For regularised deconvolution, we plot NMSE versus λ Regularisation strength0.760.780.800.820.840.860.88 N S P Extended Source: Gradient Regularisation Regularisation strength0.600.650.700.750.800.85 N S P Extended Source: Gradient Regularisation Regularisation strength0.30.40.50.60.70.8 N S P Extended Source: Gradient Regularisation
FIG. 2: NSP (of clean map generated using gradient regularization) versus λ for extended source with strong, weak and very weak signal strengths respectively. The plots show that, as expected, optimal strength of regularisation, which maximises NSP,reduces with the strength of the source one is probing. However, the curves are reasonably ‘flat’, indicating that a broad rangeof values of λ can yield near-optimal result. Source Map
Source Map
Source Map
Dirty Map -29.8804 30.8062
Dirty Map -7.92312 8.86577
Dirty Map -4.41304 4.53962
Clean Map: Unregularised, -0.0233686 0.0727715
Clean Map: Unregularised, -0.0110689 0.0237456
Clean Map: Unregularised, -0.00891336 0.0110235
Clean Map: Gradient regularisation with strength = 5 -0.0114811 0.0733919
Clean Map: Gradient regularisation with strength = 50 -0.00122189 0.0160646
Clean Map: Gradient regularisation with strength = 500
FIG. 3: Illustration of gradient regularised deconvolution for extended sources. Each column corresponds to a progressivelydecreasing strength of injected signal with a skymap that resembles the Milky Way galaxy. The rows correspond to injectedmap, dirty map, unregularised and regularised clean maps. The number of iterations for the unregularised clean maps are12 , , λ was chosen as 5 , ,
500 respectively with 100 iterations. For a strongsource, one can see that deconvolution without regularization produces good reconstruction of the true sky, regularization isless effective here. However, when the source is weak, regularisation brings dramatic improvement in deconvolution. For weaksources, the unregularised clean maps (even the dirty maps) are dominated by noise, while regularised deconvolution reducesnoise and brings some of the source map features above the noise level. Quantitative measures are presented in Table I. Regularisation strength1.01.52.02.53.03.5 N M S E Point Source: Norm Regularisation Regularisation strength010203040506070 N M S E Point Source: Norm Regularisation Regularisation strength050100150200250300 N M S E Point Source: Norm Regularisation
FIG. 4: NMSE (of clean map generated using norm regularization) versus λ for point sources with strong, weak and very weaksignal strengths respectively. The plots show that the optimal strength of regularisation that minimises NMSE, reduces withthe strength of the source. Notice that the curves are nearly ‘flat’ beyond a certain λ , indicating that a broad ranges of valuesof λ can yield near-optimal results. Source Map
Source Map
Source Map
Dirty Map -14.1178 29.838
Dirty Map -4.05032 6.3681
Dirty Map -3.57506 3.66798
Clean Map: Unregularised, -0.0136952 0.0455763
Clean Map: Unregularised, -0.0087932 0.0103909
Clean Map: Unregularised, -0.00899822 0.00935227
Clean Map: Norm regularisation with strength 100 -0.0120374 0.0418591
Clean Map: Norm regularisation with strength 1000 -0.00187679 0.00302768
Clean Map: Norm regularisation with strength 1000 -0.0018231 0.00161395
FIG. 5: Illustration of norm regularised deconvolution for a localised source. Each column corresponds to a progressivelydecreasing strength of injected sky-map with a point source at an arbitrary pixel. The rows correspond respectively to injectedmap, dirty map and clean map without and with regularisation. The number of iterations for unregularised and regulariseddeconvolution are chosen as 3 and 50 respectively. For regularised deconvolution λ was chosen as 100 , , in Fig. 4. The NMSE curve remains nearly constant for alarge range of λ ; this indicates that any choice of λ abovea certain critical value would produce nearly optimal re-sults. We set λ = 100 for the strong sources and 1000 forthe weaker sources. Regularisation improves the qualityof deconvolution for weak sources, as seen in the last rowof Fig. 5. NMSE for a very weak point source improvesby a huge factor, from 9 . . σ and 3 σ cutoff levels where σ is the standard deviation of the re-spective maps. Masking is a procedure where pixel valuesof the map are set to zero if they are smaller than thecutoff value. The masking was performed for the cleanmap corresponding to the third column of Fig. 5. Wesee that a 3 σ mask correctly localizes the injected sourcein the regularised map. The masked unregularised cleanmap picks up a wrong pixel as shown in Fig. 6. This exer-cise is useful for detecting ‘outliers’—significantly brightpixels in the skymap. If the deconvolution itself createsoutliers, the usefulness of the exercise will be severelycompromised. These results also signify that, althoughwe are using NMSE as one of the possible quality mea-sures for reconstruction (which is certainly better thanNSP for point sources), it may not capture the full ex-tent of advantage of regularisation. C. Robustness of regularised deconvolution
So far we have considered specific injections in thisstudy. We now show how any reasonable choice of theregularisation constant and number of iterations generi-cally improves the quality of deconvolution.
1. Insensitive to the number of iterations
A considerable advantage of regularisation comes inthe form of numerical stability of deconvolution. Qualityof reconstruction of the SGWB map using unregulariseddeconvolution deteriorates after a few iterations due tothe accumulation of numerical noise. This can be seen inthe top panels of Fig. 7. In the figure, NSP and NMSEare plotted against the number of iterations for a strongextended source (left) and a strong point-source (right)[same injections were used in the left column of Fig. 3 and5 respectively]. The plots for regularised deconvolutionare provided in the lower panels of Fig. 7, which showthat the quality of deconvolution stabilises after ∼ −
20 iterations for these cases. In accordance with thisresult, we choose to use 50 −
100 iterations for regulariseddeconvolution to be on the safer side (though a smallernumber could have yielded similar results).
2. Nearly insensitive to the choice of λ Although we pick a value for the regularisation con-stant λ for each injection separately such that regular-isation produces optimal results, as seen in the Fig. 2and 4, we demonstrate that a broad range of values for λ could produce similar results. Our results do not re-quire a fine-tuned value of λ . We find that any choiceof λ in a given range would produce fairly similar re-sults as long as the strength of the source is similar, i.e.,the procedure is not very sensitive to the exact shapeand features of the source intensity distribution in thesky. For instance, λ = 1000 or 10000 would produceacceptable results (within ∼
10% of the best reconstruc-tion quality) for all the six injections considered in Fig. 3and 5. It is worth noting that stronger the regularisa-tion strength, larger the bias. Also, strong gradient reg-ularisation washes out the finer details, while a strongnorm regularisation reduces the strength of the source.Although a visual inspection of our result confirms thisaspect, it may not always be reflected in the NSP andNMSE measures. Therefore in order to capture the max-imum information from the data, it is recommended thatthe smallest value of λ that produces a reasonable recon-structed map be used.
3. Simulations
We now demonstrate that the results are indeed in-sensitive to the choice of regularisation strength by per-forming simulations with a fixed choice of λ . We im-plement our regularized deconvolution routine on 1000simulations for each of the two cases - point sources andextended sources, corresponding to weak to very weaksignal strength. For each simulation, we generate a dirtymap by convolving the injected map with the beam andadd a noise map. The noise map is generated by process-ing Gaussian noise in frequency domain correspondingto the two LIGO detectors. We follow the procedure de-scribed in Mitra et al. [17].For point source injections, we randomly select a pixeland assign an intensity drawn from a uniform distribu-tion in the ‘weak’ to ‘very weak’ range. We deconvolvethe dirty map without regularisation and with norm reg-ularisation with λ = 1000. We then make histograms ofNMSE and its difference obtained from these two typesof clean maps, as presented in the right panel of Fig. 8.The process is more elaborate for extended sources.We compute the angular power spectrum, C l , of theMilky-Way-Galaxy-like map that we have been using inthis paper. We then generate simulated maps from this C l using HEALPix tools. We take the absolute valueof the maps to make all the pixels positive and finallymultiply with random scaling factors uniformly drawnfrom a suitable range; the range is chosen such that thestrength of the injections lies in the range ‘weak’ to ‘veryweak’. We then create the dirty maps and deconvolve0 Dirty Map: 2-sigma mask
Unregularised Clean Map: 2-sigma mask
Norm regularised Clean Map: 2-sigma mask
Dirty Map: 3-sigma mask
Unregularised Clean Map: 3-sigma mask
Norm regularised Clean Map: 3-sigma mask
FIG. 6: Plots of positive 2 σ (top row) and 3 σ (bottom row) outliers for dirty map (left) and clean maps with (right) andwithout (middle) regularisation, corresponding to the very weak point source injection shown in Fig. 5. Spurious localisedsources appear in all the maps with a 2 σ mask. Even with a 3 σ mask dirty map and unregularalised clean map show spuriousoutliers, while the norm regularised clean map precisely locates the source pixel. N S P Extended Source: No regularisation N M S E Point Source: No regularisation N S P Extended Source: Gradient Regularisation with strength 5 N M S E Point Source: Norm Regularisation with strength 100
FIG. 7: Stability of regularization with iteration number. The top two plots show how increasing the number of iterationsdeteriorates the quality of deconvolution for strong extended (left) and point (right) sources [these injections correspond to theleft columns of Figures 3 and 5 respectively]. Regularisation stabilises the quality of deconvolution after 10 −
20 iterations, asshown in the bottom plots for the corresponding sources with gradient and norm regularisations respectively. them without regularisation and with gradient regular-isation with λ = 500. We make histograms of NSP forthese two types of clean maps and its difference, as shownin the left panel of Fig. 8.The histograms in Fig. 8 clearly show that regu- larised deconvolution (slanted hatched bars) significantlyoutperforms unregularised deconvolution (horizontallyhatched bars). The histograms of differences (plain bars)of NSP or NMSE between with and without regularisa-tion (sign chosen appropriately) do not show any negative1values. This further unveils that regularisation improvesreconstruction not only statistically, but for each individ-ual simulation. D. Prescription for Real Data
We have now demonstrated that regularisation im-proves the quality of deconvolution and it is robustagainst the choice of parameters, provided we have someidea about the strength of the signal we are trying toprobe. Since no SGWB has been detected so far, we ex-pect that very weak sources will be the primary targetsof a search in the coming years. On the contrary, wewere unable to faithfully reconstruct skymaps that aremuch weaker than the ones considered here [third col-umn of Fig. 3 and 5] irrespective of the regularisationmethod and parameters. Given this limitation, a prac-tical approach would be to target the weakest case thatfall in the realm of validity of our procedure. From thisperspective, to apply the procedure on the data from theupcoming observing runs of the current detectors, we rec-ommend finding the optimal strength of regularisation forthe weakest sources through a set of simulated injectionsas prescribed in this work. If the normalisation of data,beams, and noise properties are similar to those used inour simulation, we expect the optimal strength to lie inthe range λ ∼ − λ over a large range of valuesand see if anything interesting stands out in the map.However, this may not be the best practice to follow forall-sky all-frequency searches [22, 24], where the numberof maps is already very high (equal to the number offrequency bins), making it computationally challengingand significantly boosting the risk of discovering ‘lookelsewhere’ effects.One could also argue if the dirty maps could be directlyused without performing deconvolution at all. In fact, forcertain cases in Table I, values of NSP or NMSE are in-deed better for dirty maps. However, there are difficultieswhich prevent us from taking such advantages. First, thequality of dirty maps is worse than the regularised mapsfor very weak sources, which are going to be our pri-mary targets. Reading from Table I, for the very weakextended source, NSP dirty = 0 . < NSP grad = 0 .
79 andfor the very weak point source, NMSE dirty = 14 . > NMSE norm = 1 .
29. Second, one needs to have a sourcesky model to use the dirty maps directly. Regulariseddeconvolution does not need a specific source model, ituses generic features of the source. In order to look fora point source in the dirty map whose location is notknown, perhaps the most reasonable way would be tofind outliers in the map, like in Fig. 6, where the regu-larised clean map outperforms the dirty map. Similarly,to find whether the dirty map in Fig. 3 has an embeddedsource map, assuming that we had a model for the map(the Milky-Way-like pattern in this case, perhaps fromelectromagnetic surveys), we could find the NSP and test its statistical significance (similar to ‘Matched Filtering’the sky). Even in such cases, one would have to worryabout the accuracy and completeness of the model. Inthe absence of a specific source model, there is no obviousway to check if the dirty map embeds an extended sourceand infer its shape, and therefore one has to work withthe clean map.
VI. CONCLUSION
Mapping an anisotropic stochastic gravitational wavebackground using data from ground-based detectors isbecoming progressively important as detectors are break-ing sensitivity barriers and new cosmological results arebeing published. The task however is challenging. Onefundamental hurdle is that the matrix that connects thesource sky to the data is somewhat ill-conditioned, mak-ing it non-trivial to deconvolve the filtered cross-spectraldata from pairs of detectors, a.k.a. the dirty map. Inthis work, we demonstrated that regularized deconvolu-tion provides a robust yet straightforward way to addressthis issue and the method can be readily applied to thecurrent LIGO-Virgo analyses.Motivated by an earlier work on gravitational lens-ing [18], we introduced and applied regularised decon-volution in SGWB mapmaking. We use two forms ofregularisation function here: (1) norm regularisation thattries to minimise power in the whole map, which is suit-able for localised sources (2) gradient regularisation thattries to reject small angular scale variations, which issuitable for extended sources. We use a Bayesian analy-sis to determine the optimal strength of regularization forthe above two functions. The merits of these regularisa-tion schemes are demonstrated using multiple examplesof different source distributions in this paper. We showthat regularisation dramatically improves the quality ofreconstruction for weak sources, which are likely to be theprimary candidates for the first detection. The method isnot sensitive to a specific choice of the strength of regu-larisation, as long as the strength is chosen from a broadrange of values determined by following the prescriptiongiven here. For strong sources regularisation is ineffec-tive and sometimes worsens the quality of reconstructioncompared to unregularised deconvolution. However, ir-respective of the strength of the signal, regularisationstabilises the quality of deconvolution against iterations,making it a safer practice to follow.In the detection problem, one is often interested toknow the presence or absence of a particular known pat-tern in the data. In such a case, the likelihood ratiois the optimal detection statistics [26], where it is as-sumed that the pattern in the sky is accurately known.However, most often in the real scenarios, only a pieceof information about the source intensity distribution isknown, and it becomes vital to be able to regulate thestrength of the prior. Implementing regularization in aBayesian framework provides this versatility.2
Histogram for extended sources: Regularisation strength 500
Gradient regularisationNo regularisationWith - Without Regularisation 2 4 6 8 10 12NMSE02004006008001000
Histogram for point sources: Regularisation strength 1000
Norm regularisationNo regularisationWithout - With Regularisation
FIG. 8: The above histograms show that regularized deconvolution gives better reconstruction compared to unregularizeddeconvolution even when the regularization constant is not fine-tuned for a specific strength or pattern. A simulation wasperformed with randomly chosen skymaps (left) or point sources injected in a random direction (right) with randomly chosenstrengths of the injection in both the cases. Gradient regularisation was performed for extended sources and norm regularisationfor point sources with fixed values of regularisation constant, λ = 500 and 1000 respectively. Regularised clean maps (slantedhatched bars) have a significantly better quality of deconvolution compared to unregularised maps (horizontally hatched bars).In fact, the differences (unhatched bars) between NSP & NMSE of regularised & unregularised maps are always positive,implying that regularisation improved the quality of deconvolution in each of the simulations, not just statistically! Our method is best suited for blind searches wherethe source location or intensity distribution is unknown,but we can use non-informative priors like gradient andnorm regularization. It is worth mentioning here thatunlike the detection problem where one is interested toknow if a particular pattern is present in the data, themapmaking problem concerns with reconstruction of themost probable pattern in the sky given the data. Onecan, in principle, use this reconstructed map to test thepresence or absence of a specific pattern, but it may bemore optimal to tailor an informative prior in such cases.The clean maps generated by our method are readily use-ful to search for anomalies and outliers caused by un-known sources or instrumental artifacts. It is particu-larly relevant for conducting blind all-sky narrowbandsearches [24], though we have not yet tested our methodfor that application, which may pose a computationalchallenge if not implemented carefully.Here we have considered only the LIGO Hanford andLivingston detectors. The method can be easily extendedto a network of detectors by replacing the beam matrixand dirty map by their sums from the individual base-lines (with carefully chosen normalisation) [17, 26, 27].The network acts as a natural regulariser to a certainextent [26, 27], so it may require a smaller regularisationstrength for the same pixel resolution. However, whenKAGRA and LIGO-India join the network, the base-line lengths will significantly increase, leading to a muchhigher resolving power. This will demand finer pixels (in-creasing the number of pixels by a factor of 4 or more),making it challenging to invert the beam matrix, possiblynecessitating an even stronger regularisation.The stochastic searches are routinely being conductedin Spherical Harmonic basis [28]; deconvolution is a big challenge in this basis as well [27]. While for the caseof pixel-based radiometer analysis no regularised decon-volution has been conventionally used, for the case ofspherical harmonic analysis a heuristic SVD-based regu-larisation has been implemented. Based on the advantagewe gain by using Bayesian regularization in pixel basis,we are encouraged to propose a similar approach to thespherical harmonic basis as well.In the recent past, many theoretical models havebeen introduced to predict anisotropic SGWB from com-pact binaries [29, 29–32] and Nambu-Goto CosmicStrings [33]. These models predict an angular powerspectrum of the sky, C l , with part of the power accessiblein the sensitive frequency band of the ground-based laserinterferometric detectors, which will allow constraining ofthese models using data from those detectors through theestimation of C l [28]. Since the dirty maps are convolu-tion of the source with an extended asymmetric beam, es-timation of C l of the true sky, will require a reliable cleanmap where our method should prove to be handy. Theimprovement of NSP through regularisation for extendedsources should help studies performing cross-correlationbetween SGWB maps and the large scale structures. Go-ing beyond, since the regularization scheme applied inthis paper for SGWB mapmaking is fairly generic, wehope that it will open new avenues to aid the above men-tioned and other related investigations in the future. Acknowledgments
We would like to thank Sherry H. Suyu and VukMandic for suggesting and discussing this possibilitymany years ago! We also thank Joe Romano for his valu-3able comments. This work significantly benefitted fromthe interactions with the Stochastic Working Group ofthe LIGO-Virgo Scientific Collaboration. We acknowl-edge the use of IUCAA LDAS cluster Sarathi for thecomputational/numerical work. This research benefitedfrom a grant awarded to IUCAA by the Navajbai RatanTata Trust (NRTT). SM acknowledges support from theDepartment of Science and Technology (DST), India pro-vided under the Swarna Jayanti Fellowships scheme. JSacknowledges the support by JSPS KAKENHI GrantNumber JP17H06361. This article has a LIGO documentnumber LIGO-P1800278 and IUCAA Preprint numberIUCAA-07/2019. SB acknowledges the financial sup-port provided by the SUGWG group at the SyracuseUniversity under the grant number PHY-1707954 andthe DARKGRA group in La Sapienza under the Eu-ropean Union’s H2020 ERC, Starting Grant agreementno. DarkGRA–757480 and support from the Amaldi Re-search Center funded by the MIUR program “Diparti-mento di Eccellenza” (CUP: B81I18001170001).
Appendix A: Derivation of expression for S mp The most probable solution is defined as, ∇ S M ( S mp ) = ∇ S χ ( S mp ) + µ ∇ S R ( S mp ) = 0 . (A1)In this section we prove that the S mp can be calculatedas S mp = ( B + λ C ) − D .Let R ( S ) be any quadratic function of S with a min-ima located at S reg and let C be the hessian. By Taylorexpanding about the minima S reg , R ( S ) can be writtenas, R ( S ) = S reg + 12 S T CS (A2)For both the regularization functions we consider in ourstudy, S reg = 0, and, therefore, R ( S ) = 12 S T CS . (A3)Now, recall that N = κ B and λ = µκ . Then, χ ≡
12 (
BS − D ) T N − ( BS − D )= µ λ ( BS − D ) T B − ( BS − D ) . (A4) Next, we compute the functional form of ∇ S M ( S ). ∇ S M ( S ) = ∇ S χ ( S ) + µ ∇ S R ( S )= µ λ ∇ S (cid:2) ( BS − D ) T B − ( BS − D ) (cid:3) + µ ∇ S (cid:2) S T CS (cid:3) = ( µ/λ )( BS − D ) + µ ( CS ) . (A5)To solve for S mp , we set ∇ S M ( S mp ) = 0 and get, BS mp + λ C S mp = D , (A6)that is, S mp = ( B + λ C ) − D , (A7)the most probable solution. Appendix B: Bias Introduced by Regularization
The map reconstructed using regularization is a biasedestimator while the un-regularized reconstruction is un-biased. Bias ℘ is defined, ℘ = S true − < ˆ S > . (B1)For unregularized reconstruction, < ˆ S > = ( B T N − B ) − ( B T N − B ) S true = S true . (B2)From the above equation, it is clear that ℘ = 0 for theunregularized clean map. Regularization introduce thebias in the reconstructed map by altering the pixel-to-pixel covariance matrix of the reconstructed map.For the case of regularized deconvolution, fromEq. (A7) the bias ℘ can be calculated as, < S mp > = ( B + λ C ) − < B·S + n > = ( B + λ C ) − B·S true , (B3)Hence, the bias ℘ for this case is, ℘ = S true − < ˆ S > = ( − ( B + λ C ) − B ) · S true . (B4)In the case where the source intensity is strong, regu-larization unnecessarily introduces a bias, and therefore,an unregularized deconvolution performs slightly better.The bias cannot be corrected, because it requires theknowledge of the true sky which is being estimated. [1] B. Allen, in Relativistic Gravitation and GravitationalRadiation , edited by J.-A. Marck and J.-P. Lasota (1997),p. 373, gr-qc/9604033.[2] B. Allen and A. C. Ottewill, Phys. Rev. D , 545(1997), gr-qc/9607068.[3] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Ad-hikari, V. B. Adya, et al., Physical Review Letters ,091101 (2018), 1710.05837.[4] S. Dhurandhar, H. Tagoshi, Y. Okada, N. Kanda,and H. Takahashi, Phys. Rev. D , 083007 (2011),1105.5842. [5] J. D. Romano and N. J. Cornish, Living Reviews in Rel-ativity , 2 (2017), ISSN 1433-8351.[6] S. A. Hughes, Physics of the Dark Universe , 86(2014), ISSN 2212-6864, dARK TAUP2013, URL .[7] G. Cusin, R. Durrer, and P. G. Ferreira, ArXiv e-prints(2018), 1807.10620.[8] N. Mazumder, S. Mitra, and S. Dhurandhar, Phys. Rev.D , 084076 (2014), 1401.5898.[9] P. F. Michelson, , 933 (1987).[10] N. Christensen, Phys. Rev. D46 , 5250 (1992).[11] E. E. Flanagan, Phys. Rev.
D48 , 2389 (1993).[12] B. Allen and A. C. Ottewill, Phys. Rev.
D56 , 545 (1997).[13] B. Allen and J. D. Romano, Phys. Rev.
D59 , 102001(1999), gr-qc/9710117.[14] A. Lazzarini and J. Romano, Internal working noteLIGO-T040089-00-Z, Laser Interferometer GravitationalWave Observatory (LIGO) (2004).[15] A. Lazzarini and R. Weiss, Internal working note LIGO-T040140-00-Z, Laser Interferometer Gravitational WaveObservatory (LIGO) (2004).[16] S. W. Ballmer, Class. Quant. Grav. , S179 (2006), gr-qc/0510096.[17] S. Mitra, S. Dhurandhar, T. Souradeep, A. Lazzarini,V. Mandic, et al., Phys.Rev. D77 , 042002 (2008),0708.2728.[18] S. H. Suyu, P. J. Marshall, M. P. Hobson, and R. D.Blandford, , 983 (2006).[19] M. Rakhmanov, Class. Quant. Grav. , S673 (2006),gr-qc/0604005. [20] A. C. Searle, P. J. Sutton, M. Tinto, and G. Woan, Class.Quant. Grav. , 114038 (2008), 0712.0196.[21] M. Ryle and A. Hewish, , 220 (1960).[22] A. Ain, J. Suresh, and S. Mitra, Phys. Rev. D98 , 024001(2018), 1803.08285.[23] A. Ain, P. Dalvi, and S. Mitra, Phys. Rev.
D92 , 022003(2015), 1504.01714.[24] E. Thrane, S. Mitra, N. Christensen, V. Mandic, andA. Ain, Phys. Rev. D , 124012 (2015).[25] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt,F. K. Hansen, M. Reinecke, and M. Bartelman, Astro-phys. J. , 759 (2005), astro-ph/0409513.[26] D. Talukder, S. Mitra, and S. Bose, Phys. Rev. D83 ,063002 (2011), 1012.4530.[27] E. Thrane, S. Ballmer, J. D. Romano, S. Mitra,D. Talukder, S. Bose, and V. Mandic, Phys. Rev. D ,122002 (2009).[28] B. P. Abbott et al. (LIGO Scientific, Virgo) (2019),1903.08844.[29] A. C. Jenkins, M. Sakellariadou, T. Regimbau, andE. Slezak, Phys. Rev. D98 , 063501 (2018), 1806.01718.[30] G. Cusin, I. Dvorkin, C. Pitrou, and J.-P. Uzan, Phys.Rev. Lett. , 231101 (2018), 1803.03236.[31] G. Cusin, C. Pitrou, and J.-P. Uzan, Phys. Rev.
D97 ,123527 (2018), 1711.11345.[32] G. Cusin, C. Pitrou, and J.-P. Uzan, Phys. Rev.
D96 ,103019 (2017), 1704.06184.[33] A. C. Jenkins and M. Sakellariadou, Phys. Rev.