Anisotropies in the astrophysical gravitational-wave background: Predictions for the detection of compact binaries by LIGO and Virgo
Alexander C. Jenkins, Mairi Sakellariadou, Tania Regimbau, Eric Slezak
KKCL-PH-TH/2018-23
Anisotropies in the astrophysical gravitational-wave background:Predictions for the detection of compact binaries by LIGO and Virgo
Alexander C. Jenkins ∗ and Mairi Sakellariadou † Theoretical Particle Physics and Cosmology Group, Physics Department,King’s College London, University of London, Strand, London WC2R 2LS, United Kingdom
Tania Regimbau ‡ Universit´e Savoie Mont Blanc, CNRS/IN2P3, Laboratoire d’Annecy-le-Vieuxde Physique des Particules (LAPP), 74941 Annecy, France andUniversit´e Cˆote d’Azur, Observatoire de la Cˆote d’Azur, CNRS,Laboratoire Artemis, CS 34229, 06304 Nice Cedex 4, France
Eric Slezak § Universit´e Cˆote d’Azur, Observatoire de la Cˆote d’Azur, CNRS,Laboratoire Lagrange, CS 34229, 06304 Nice Cedex 4, France (Dated: 11th September 2018)We develop a detailed anisotropic model for the astrophysical gravitational-wave background,including binary mergers of two stellar-mass black holes, two neutron stars, or one of each, whichare expected to be the strongest contributions in the LIGO-Virgo frequency band. The angularspectrum of the anisotropies, quantified by the C (cid:96) components, is calculated using two complementaryapproaches: (i) a simple, closed-form analytical expression, and (ii) a detailed numerical study usingan all-sky mock light cone galaxy catalogue from the Millennium simulation. The two approachesare in excellent agreement at large angular scales, and differ by a factor of order unity at smallerscales. These anisotropies are considerably larger in amplitude than e.g. those in the temperature ofthe cosmic microwave background, confirming that it is important to model these anisotropies, andindicating that this is a promising avenue for future theoretical and observational work. I. INTRODUCTION
The first direct detections of gravitational waves (GW)by the Advanced LIGO [1, 2] and Advanced Virgo [3]detectors, from the inspiral and merger of pairs of blackholes [4–7], and the recent first detection of GWs fromthe inspiral and merger of a pair of neutron stars [8],have opened a new window to the Universe, its physicalprocesses, and the astrophysical sources that populate it.Each of these detections are associated with individualloud events, but it is expected that there are many morequiet compact binary mergers in the Universe that are toofaint or too distant to be individually resolved. Signalsfrom these quiet events superimpose to create a stochasticgravitational-wave background (SGWB), which is discern-ible from instrumental noise by cross-correlating the out-put from multiple detectors, and which is expected to bedetected by Advanced LIGO and Advanced Virgo aftera few years of operation at design sensitivity [4]. Manyother sources, both astrophysical (e.g., rotating neutronstars, supernovae, core collapse formation of black holesor neutron stars) and cosmological (e.g., phase transitions,inflation, cosmic strings), are expected to contribute tothe SGWB. The cosmological background can provide ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] unique information about the early stages of our Universe,while the astrophysical background probes the Universe’sevolution since the beginning of star formation.Gravitational-wave sources with an inhomogeneous spa-tial distribution lead to a SGWB characterized by pre-ferred directions, and therefore anisotropies. This is ana-logous to the measured temperature anisotropies in thecosmic microwave background (CMB) radiation. Methodsto detect these SGWB anisotropies using radiometer andspherical harmonic techniques have been developed [9–16]and applied to data from Advanced LIGO’s first observa-tional run to set upper limits on the energy flux and energydensity of the SGWB as a function of sky position [17].Recently, Refs. [18–20] have investigated a general frame-work to model the anisotropies in the SGWB induced byvarious types of sources, which was applied to the case ofcosmic strings in Ref. [20] and to the case of binary blackhole (BBH) mergers in Ref. [21].In what follows we use a result [Eq. (26)] of Ref. [20]to study anisotropies induced in the SGWB by astrophys-ical sources, namely from compact object mergers withinanisotropically distributed galaxies. We assume that thebackground is stationary, which implies that the statist-ical properties of the SGWB do not vary over the timescale of the observation. Our analysis follows two distinctapproaches. In the former (analytical approach), we usesimple analytical functions for the galaxy number densityand the galaxy-galaxy two-point correlation function toderive a simple, closed-form expression for the two-pointcorrelation function and its multipole components C (cid:96) . Let a r X i v : . [ a s t r o - ph . C O ] S e p us point out that while Ref. [21] also used a two-pointcorrelation function, our study differs in several ways: (i)we include all three different types of merger: binariesconsisting of two neutron stars, two stellar-mass blackholes, or one of each; (ii) we calculate the kinematic di-pole, which we then subtract since it interferes with theanisotropy statistics (as demonstrated in Ref. [20]); (iii)we follow the fiducial model of the LIGO/Virgo collab-oration paper [22] rather than the astrophysical modelof Refs. [23–25] to estimate the rate of mergers; (iv) weuse a nonlinear power law expression Eq. (64) to modelthe two-point galaxy clustering, rather than basing thison linear transfer functions of the matter overdensity asin Ref. [21]. In the latter (catalogue approach), we makeuse—for the first time in the SGWB literature—of a real-istic mock catalogue of galaxies and employ recipes toinfer the production and merger rates of compact objects.In both approaches, we restrict our attention to themost important sources for LIGO, Virgo, and otherplanned ground-based detectors. Future space-based de-tectors such as the Laser Interferometer Space Antenna(LISA) [26] will probe a lower-frequency window of theGW spectrum, and will therefore be sensitive to differentpopulations of astrophysical systems, which will producea different stochastic background signal. However, muchof Sec. II is still valid in this case, and we intend to repeatour analysis for the LISA frequency band in a future work.Our present study is organized as follows. In Sec. IIwe derive a general expression for the astrophysical grav-itational wave background (AGWB) induced by an an-isotropic distribution of galaxies hosting compact binarymergers. We first write the SGWB density parameter Ω gw in terms of the average number density of galaxiesper comoving volume, the galaxy number overdensity, therate of binary mergers per galaxy, the observer’s peculiarvelocity, and the GW strain spectrum of each binary mer-ger. These quantities depend on the star formation rates(SFRs) and metallicities of the galaxies, and the spinsand masses of the two binary components. We performour analysis using the fiducial model of Ref. [22], andthus consider the merger rate, the distribution of binaryparameters, and the gravitational waveforms emitted bythe binary mergers as those of Ref. [22]. We then decom-pose the density parameter as the average (isotropic) part¯ Ω gw and the GW density contrast δ gw , which encodesthe anisotropies. In Sec. III we use analytical functionsto describe the average galaxy number density and thegalaxy-galaxy two-point correlation function and derivea simple closed-form expression for each of the C (cid:96) coef-ficients, up to a frequency-dependent factor which weevaluate numerically. In Sec. IV we follow a more ac-curate approach and use a realistic mock catalogue ofgalaxies, hence relaxing some of the assumptions madein the analytical approach. In Sec. V we compare theanalytical and numerical approaches, and comment on theimprint to the SGWB from the anisotropic distributionof compact object mergers as compared to those inferredby a cosmic string network [20]. Some technical details are given explicitly in the appendixes. II. GENERAL EXPRESSION FOR THE AGWB
The dimensionless density parameter expressing theintensity of a SGWB, with observed frequency between ν o and ν o + d ν o and arriving from an infinitesimal solidangle d σ o centered on the direction ˆ e o , is defined as Ω gw ( ν o , ˆ e o ) ≡ ρ c d ρ gw d(ln ν o ) d σ o = 8 π Gν o H d ρ gw d ν o d σ o , (1)where we have used the customary normalization withrespect to the critical density ρ c = 3 H / (8 π G ). Tostudy anisotropies in the SGWB induced by astrophysicalsources, we consider a Friedmann-Lemaˆıtre-Robertson-Walker spacetime, and neglect cosmological perturbations,keeping only the anisotropy due to the source density con-trast and the dipole induced by the peculiar motion ofthe observer. Following Ref. [20], we thus have Ω gw ( ν o , ˆ e o ) = π t H ν o ) (cid:90) ∞ d z zE ( z ) × (cid:90) d ζ ¯ nR (1 + δ n + ˆ e o · v o ) (cid:90) S d σ s r ˜ h , (2)where t H ≡ /H o is the Hubble time, z is the redshift, ζ represents the set of source parameters, ¯ n is the average(homogeneous) source number density per comoving unitvolume, R is the rate of gravitational wave bursts persource, δ n ≡ ( n − ¯ n ) / ¯ n is the source number overdensity, ˆ e o is the observation direction, v o is the observer’s pe-culiar velocity, ˜ h is the GW strain spectrum of a burst,and the final integral is over a sphere centered on thesource. Note that Eq. (2) has been modified with respectto Eq. (26) of Ref. [20] by integrating over redshift ratherthan conformal time, and by using the number densityper comoving volume, rather than per physical volume.Assuming the standard flat Λ CDM cosmology, we have E ( z ) ≡ H ( z ) H o = (cid:113) Ω m (1 + z ) + Ω Λ , (3)with Ω m = 0 . Ω Λ = 0 . H o =67 . − Mpc − .Note that, while Eq. (26) of Ref. [20] imposed a cutofftime on the line-of-sight integral to remove nearby, loud,infrequent, resolvable sources and isolate the stochasticpart of the signal (see Sec. IIC of Ref. [20]), this is not donein Eq. (2) above. This is because the background frombinary mergers consists of far fewer GW-emitting eventsthan the cosmic string background considered in Ref. [20].As a result, the problem of distinguishing different eventsfrom each other is greatly reduced, and the detector noisebecomes the main barrier to resolving events individually.Indeed, a network of future “third-generation” detectorsis expected to be able to resolve the vast majority ofbinary mergers in this frequency window, including morethan 99.9% of BBH mergers [22, 27, 28]. (For a discussionof how such a large ensemble of individually resolvedmergers may be used as a cosmological probe, see Refs. [29,30].) The astrophysical background from unresolvablebinary mergers therefore depends on the sensitivity ofthe detector network, and future detectors will be able togreatly reduce the level of the background by resolvingbinaries out to much larger redshifts. In order to phraseour results in a detector-independent way, we have optedto include all binary merger events, no matter how closethey are to the observer.In what follows, we study the imprint of an anisotropicdistribution of binary mergers, within galaxies, on theSGWB. Since any anisotropies on scales smaller than thetypical size of a galaxy are inaccessible to us, we treatgalaxies as point sources. We therefore interpret n as thenumber density of galaxies, and R as the rate of binarymergers per galaxy. Let us write ζ = (cid:0) ζ b , ζ g (cid:1) , (4)where ζ g are the parameters of the galaxy (mass, age,luminosity, metallicity, etc.) and ζ b are the parameters ofthe compact binary (masses and spins of the components).We then have n = n (cid:0) z, ˆ e o , ζ g (cid:1) , ˜ h = ˜ h ( ν s , ˆ e s , ζ b ) , R = R (cid:0) z, ζ g , ζ b (cid:1) , (5)since the galaxy number density is independent of thesource parameters, the strain is independent of the galaxyparameters, and the rate per galaxy depends on both.Note that ˜ h is a function of the source-frame frequency ν s , which is related to the observed frequency ν o by ν s = ν o (1 + z )[1 + ˆ e o · ( v g − v o )] , (6)where v g and v o are the peculiar velocities of the galaxyand the observer, respectively. We allow each source toemit GWs anisotropically, so that ˜ h is also a functionof the direction of propagation away from the source’sposition, ˆ e s .We can capture much of the relevant astrophysicalinformation with just six parameters: the SFR ψ andmetallicity Z of the galaxies, and the masses and spinsof the two binary components. The compact objects weconsider are all the end products of stellar evolution, sowe can use the SFR history of a galaxy to calculate itspopulation of compact binaries, and therefore parameter-ize the rate of mergers of these objects. The formation ofmassive black holes from stars is inhibited by stellar windsin high-metallicity environments, so we also include thegalaxy metallicity as a parameter to account for this. Allother parameters describing the galactic environment areneglected. We further assume that the compact binaryorbits are circular and have nonprecessing spins (this isexpected to be the case for mergers resulting from isolatedbinary evolution—such binaries are expected to circularizelong before merger due to gravitational backreaction). Let ψ ( z ) be the SFR of a galaxy at redshift z in units of M (cid:12) yr − , and Z be the metallicity of the galaxy (i.e. thefraction of the galaxy’s mass that is in elements heavierthan helium). In order to use the SFR to parameterizethe merger rate, we must take into account the time delaybetween stars being formed, evolving to become compactobjects, and eventually merging with each other. Thesedelay times are typically much longer than the timescalesover which the SFR of a galaxy changes, so cannot beneglected. We therefore define the delayed SFR, ψ d ( z ) ≡ (cid:90) ∞ d t d p ( t d ) ψ ( z f ) , (7)where z f ( z, t d ) is the redshift at which the stars are formed,at a time t d before the merger occurs at redshift z . Thisis the convolution of the SFR with the probability dis-tribution for the delay times, which varies depending onthe objects we consider. Certainly, not all stars becomemerging compact objects, so ψ d is not equal to the mergerrate—but it is proportional to it. It is therefore muchsimpler to use ψ d as a galaxy parameter, rather than ψ .Similarly, rather than using the metallicity Z , it is moreconvenient to use a logarithmic scaling, and to normalizerelative to the solar metallicity Z (cid:12) ≈ .
02. So we define
Z ≡ log ZZ (cid:12) , Z ∈ ( −∞ , Z max ] , (8)where Z max ≡ − log Z (cid:12) ≈ . ζ g = ( ψ d , Z ) , ζ b = ( m , m , χ , χ ) , (9)where ψ d and Z are the delayed SFR and log-normalizedmetallicity introduced above, m and m are the massesof the two compact objects in a binary (measured in unitsof M (cid:12) ), and χ ≡ S m , χ ≡ S m , χ , χ ∈ [ − , +1] , (10)are the dimensionless spin parameters of the compactobjects (where S i is the spin angular momentum). Usingthese parameters, we consider three different types ofmerger: binaries consisting of two neutron stars, twostellar-mass black holes, or one of each. We call theseBNS, BBH, and BHNS respectively. Each type of binaryhas a different average merger rate, a different distributionover the parameters ζ b , and a different distribution forthe delay time t d (giving a different delayed SFR ψ d ).Equation (2) therefore becomes Ω gw ( ν o , ˆ e o ) = (cid:88) i π t H ν o ) (cid:90) z max d z zE ( z ) × (cid:90) d ζ g ¯ n (1 + δ n + ˆ e o · v o ) (cid:90) d ζ b R i S i , (11)where index i runs over BBH, BNS, and BHNS, and forbrevity we have introduced the variable S i ( ν s , ζ b ) ≡ (cid:90) S d σ s r ˜ h i . (12)We have set an upper limit z max = 10 on the redshiftintegral, reflecting the fact that the merger rate is essen-tially 0 at redshifts greater than this. The galaxy andbinary parameter integration measures are (cid:90) d ζ g = (cid:90) + ∞ d ψ d ,i (cid:90) Z max −∞ d Z , (13) (cid:90) d ζ b = (cid:90) ∞ d m (cid:90) ∞ d m (cid:90) +1 − d χ (cid:90) +1 − d χ , (14)where we have allowed the delayed SFR to be differentfor different types of binary, by allowing a different distri-bution of delay times for each case, ψ d ,i ( z ) ≡ (cid:90) ∞ d t d p i ( t d ) ψ ( z f ) . (15)For the Λ CDM cosmology we consider, the formationredshift z f of an object that takes a time t d to eventuallymerge at redshift z is given by1 + z f ( z, t d )= (1 + z ) (cid:34) cosh (cid:32) Ω / Λ t d t H (cid:33) − E ( z ) Ω / Λ sinh (cid:32) Ω / Λ t d t H (cid:33)(cid:35) − / . (16) A. Intragalactic details of the model
In what follows we describe in detail how we model in-tragalactic quantities (i.e., those that are relevant withineach galaxy)—the merger rate, the distribution of binaryparameters, and the gravitational waveforms emitted bythe binary mergers. This is almost identical to the fidu-cial model in Ref. [22], to which we refer the reader for amore thorough discussion and justification of the variousmodeling choices. The only differences from Ref. [22] are(i) the inclusion of BHNS mergers, for completeness, and(ii) the use of the merger rate per galaxy, rather than percomoving volume (this simplifies our later analysis regard-ing the anisotropies in the background). The modeling ofintergalactic quantities—the galaxy number density n andits clustering statistics—is different in our two approaches,and is discussed in Secs. III and IV.First, in order to calculate S i , we use the hybrid wave-form models of Refs. [31, 32]. These are valid for BBH, and when integrated over the sphere give S BBH ≡ (cid:90) S d σ s r ˜ h = 5( G M ) / π / × ν − / (cid:104) (cid:80) i =2 α i ( π GM ν s ) i/ (cid:105) , ν s < ν c ν − / (cid:104) (cid:80) i =1 (cid:15) i ( π GM ν s ) i/ (cid:105) , ν ≤ ν s < ν c (cid:20) (cid:16) ν s − ν ν (cid:17) (cid:21) − , ν ≤ ν s < ν (17)Here we have defined the total mass and chirp mass, M ≡ m + m , M ≡ ( m m ) / M / . (18)Recall that the source-frame frequency ν s is a function ofthe redshift and the peculiar velocities [as given by Eq. (6)],which means that S i implicitly depends on these as well.The frequencies ν , ν , ν , ν are numerical constants foreach system, given by ν i = 1 π GM ν (0) i ( χ ) + (cid:88) j =1 3 − j (cid:88) k =0 y ( jk ) i (cid:18) M M (cid:19) j/ χ k , (19)where the spin parameter χ is defined as χ ≡ m M χ + m M χ . (20)For the BNS and BHNS cases, we truncate the wave-form at the merger frequency ν and ignore any higherfrequencies, since the merger and ringdown phases of theabove expression are only valid for BBH. (BNS mergersoccur at frequencies above the LIGO-Virgo band any-way, so for observational purposes this has little effect.However, it should be possible to extend our model tohigher frequencies by using waveforms that account forneutron-star (NS) matter effects. This may be desirableas future detectors improve our sensitivity to GWs atthese frequencies.) The spectrum is therefore determinedfor all types of binaries by the numerical constants α i , (cid:15) i , y ( jk ) i (which are derived from fits between numerical sim-ulations and post-Newtonian expansions), and c i (whichare chosen to ensure the waveform is continuous). Thevalues used for each of these match those in Refs. [31, 32],and are given in Appendix B.For the binary merger rate per galaxy, R i , we have R i (cid:0) z, ζ g , ζ b (cid:1) = (cid:15)f Z p i ( ζ b ) ψ d ,i , (21)where p i ( ζ b ) = p i ( m , m , χ , χ ) is the joint probabilitydistribution for the masses and spins for a compact binaryof type i . The fact that not all stars ultimately end upin compact binaries is accounted for by multiplying thedelayed SFR by the efficiency factor (cid:15) (which is someunknown constant), as well as the factor f Z , which ac-counts for the suppression of massive BH formation inhigh-metallicity environments. The factor f Z enforcesthe assumption that black holes with mass greater than30 M (cid:12) can only be formed in galaxies with Z ≤ Z (cid:12) = ⇒ Z ≤ log ≈ − . . (22)Hence one considers a rate correction factor, f Z ( Z , m , m ) = (cid:40) , m , m < M (cid:12) Θ (cid:0) log
10 12 − Z (cid:1) , else (23)where Θ ( x ) is the Heaviside step function. The delaytime distribution is modeled as p ( t d ) ∝ /t d betweensome minimum delay time t min ,i (20 Myr for BNS, 50Myr for BBH and BHNS) and the maximum t max equalto the age of the Universe at redshift z , namely t ( z ) = 2 t H Ω / Λ arcsinh (cid:32)(cid:115) Ω Λ Ω m (1 + z ) (cid:33) . (24)At redshift zero this reduces to t o ≈ . t H , the current age of the Universe. The delayed SFR is therefore ψ d ,i = 1ln ( t ( z ) /t min ,i ) (cid:90) t ( z ) t min ,i d(ln t d ) ψ ( z f ) , (25)where the prefactor gives the appropriate normalizationof the probability distribution.It remains to specify the source parameter distributions p i ( ζ b ). We take neutron stars as having masses uniformlydistributed between 1 M (cid:12) and 2 M (cid:12) , with zero spin. Thelatter is motivated by pulsar observations and by para-meter estimation on GW170817, which both indicate thatBNS systems typically have low spins ( χ (cid:46) .
05) by thetime they merge [33] (in any case, we find that spin haslittle effect on the final results). For black holes, theprimary mass is given by a Salpeter initial mass function ∝ m − . between 5 M (cid:12) and 95 M (cid:12) . If the secondary massis a black hole, then it is distributed ∝ / ( m − M (cid:12) ).We require the sum of the masses to be less than orequal to 100 M (cid:12) , and take the black hole spins as uniformbetween − p BNS ∝ δ ( χ ) δ ( χ ) ,p BBH ∝ m − . / ( m − M (cid:12) ) ,p BHNS ∝ m − . δ ( χ ) , M (cid:12) ≤ m ≤ M (cid:12) , M (cid:12) ≤ m ≤ m ≤ M (cid:12) , M (cid:12) ≤ m ≤ M (cid:12) , M (cid:12) ≤ m ≤ M (cid:12) ,m + m ≤ M (cid:12) , M (cid:12) ≤ m ≤ M (cid:12) , (26)with a normalising constant chosen appropriately in eachcase such that (cid:82) d ζ b p i ( ζ b ) = 1.The rate R i is not yet fully determined, as the con-stant (cid:15) is unknown. However, we can eliminate this con-stant using the local (i.e. redshift zero) rates inferred byLIGO/Virgo. We do this by calculating the total rate ofmergers of type i per unit comoving volume at redshift z , R i ( z ) ≡ (cid:90) d ζ g ¯ n (cid:90) d ζ b R i (cid:0) z, ζ g , ζ b (cid:1) = (cid:15) (cid:90) d ζ g ¯ n (cid:90) d ζ b f Z p i ( ζ b ) ψ d ,i . (27)We require that in the limit z → R i (0) = R (local) i . Themost up-to-date values are [22] R (local)BNS = 1 . × − Mpc − yr − , R (local)BBH = 1 . × − Mpc − yr − . (28)LIGO and Virgo have not yet observed any BHNS events,so it is only possible to place an upper limit of [34] R (local)BHNS ≤ . × − Mpc − yr − . (29) Based on previous population synthesis studies, the trueBHNS rate is expected to be somewhere between theBNS and BBH rates (e.g., see Ref. [35] for a review), andtherefore roughly an order magnitude less than this upperlimit. In order to remain agnostic about the true valueof the BHNS rate, we consider two cases: one where weinclude BHNS mergers at the maximal rate Eq. (29), andone where we set the BHNS rate to 0.Matching to these local rates allows us to eliminate (cid:15) to find R i = R (local) i I i p i ( ζ b ) f Z ψ d ,i , (30)where we have defined the normalizing constants I i = (cid:90) d ζ g (¯ nψ d ,i ) | z =0 (cid:90) d ζ b p i ( ζ b ) f Z . (31)This is particularly simple in the BNS case, since f Z = 1,so I BNS = (cid:82) d ζ g (¯ nψ d ,i ) | z =0 .Our expression for the SGWB density parameter there-fore becomes Ω gw = (cid:88) i π R (local) i I i ( t H ν o ) (cid:90) z max d z zE ( z ) × (cid:90) d ζ g ψ d ,i ¯ n (1 + δ n + ˆ e o · v o ) (cid:90) d ζ b p i ( ζ b ) f Z S i . (32) B. Decomposing the background
We can decompose the density parameter Ω gw as Ω gw ≡ ¯ Ω gw (1 + δ gw ) , (33)where ¯ Ω gw is the average (isotropic) value of the dens-ity parameter over the sky, and δ gw is the GW density contrast, which encodes the anisotropies. The latter canitself be decomposed as δ gw ≡ δ (s)gw + D ˆ e o · ˆ v o , (34)where δ (s)gw is the density contrast due to the true cosmo-logical anisotropies, and the latter term is the kinematicdipole, with direction ˆ v o ≡ v o / | v o | and magnitude givenby the dipole factor D . This factor can be calculated byperforming a Taylor expansion in x ≡ ˆ e o · v o around x = 1, giving D ( ν o ) ≡ v o ¯ Ω − ∂Ω gw ∂x (cid:12)(cid:12)(cid:12)(cid:12) x =1 ,δ n =0 , (35)where v o ≡ | v o | . For our particular case, we thereforehave¯ Ω gw = (cid:88) i π R (local) i I i ( t H ν o ) (cid:90) z max d z zE ( z ) (cid:90) d ζ g ψ d ,i ¯ n (cid:90) d ζ b p i ( ζ b ) f Z S i , (36) δ (s)gw = ¯ Ω − (cid:88) i π R (local) i I i ( t H ν o ) (cid:90) z max d z zE ( z ) (cid:90) d ζ g ψ d ,i ¯ nδ n (cid:90) d ζ b p i ( ζ b ) f Z S i , (37) D = v o ¯ Ω − (cid:88) i π R (local) i I i ( t H ν o ) (cid:90) z max d z zE ( z ) (cid:90) d ζ g ψ d ,i ¯ n (cid:90) d ζ b p i ( ζ b ) f Z ∂∂x ( x S i ) . (38)The anisotropies are then characterized by the over-density field δ (s)gw , either directly or in terms of its statistics.One particularly useful statistical descriptor is the two-point correlation function (2PCF), defined as the secondmoment of the overdensity field, C gw ( θ o , ν o ) ≡ (cid:68) δ (s)gw ( ν o , ˆ e o ) δ (s)gw ( ν o , ˆ e (cid:48) o ) (cid:69) , (39)where θ o ≡ cos − ( ˆ e o · ˆ e (cid:48) o ), and the angle brackets denotean averaging over all pairs of directions ˆ e o , ˆ e (cid:48) o whose angleof separation is θ o , as well as an ensemble averaging overpossible random realizations of the SGWB. The firstmoment (i.e., mean) vanishes by definition, and if thebackground is a Gaussian random field (GRF) then allhigher moments either vanish or are expressed in termsof the second moment by Wick’s theorem. The 2PCFtherefore uniquely characterizes the anisotropies in theGaussian part of the background. It is common prac-tice (particularly in the CMB literature) to perform amultipole expansion of the 2PCF, C gw ( θ o , ν o ) = ∞ (cid:88) (cid:96) =0 (cid:96) + 14 π C (cid:96) ( ν o ) P (cid:96) (cos θ o ) , (40)where P (cid:96) ( x ) denotes the (cid:96) th Legendre polynomial. The anisotropies are then described in terms of the C (cid:96) com-ponents, which are given by C (cid:96) ( ν o ) ≡ π (cid:90) +1 − d(cos θ o ) C gw ( θ o , ν o ) P (cid:96) (cos θ o ) . (41)The quantity (cid:96) ( (cid:96) + 1) C (cid:96) / π is roughly the contributionto the anisotropic variance of δ (s)gw per logarithmic bin in (cid:96) , as can be seen by consideringvar (cid:16) δ (s)gw (cid:17) = C gw ( θ o = 0) = (cid:88) (cid:96) (cid:96) + 14 π C (cid:96) ≈ (cid:90) d(ln (cid:96) ) (cid:96) ( (cid:96) + 1) C (cid:96) π . (42)While the overdensity δ (s)gw is a random field, the C (cid:96) ’sare treated as deterministic quantities, averaged over someensemble of realizations of the SGWB (e.g., as observedat different locations in the Universe). As a result, anymeasurement of the C (cid:96) ’s has some cosmological uncer-tainty associated with it, due to the fact that we can onlyaccess a single realization of the SGWB. This is calledthe cosmic variance, and is given byvar( C (cid:96) ) = 22 (cid:96) + 1 C (cid:96) . (43) C. Non-Gaussianity in the AGWB
As mentioned above, the analysis of the SGWB is simpli-fied if it is a GRF, as this eliminates the need for anythingother than the second moment of the overdensity (i.e. the2PCF). In Ref. [20], it was found that a GW backgroundcomposed of independent sources and discretized into N pix pixels on the sky is a GRF at frequency ν if it satisfies νT (cid:29) N pix Λ + N Λ , (44)where T is the observation time and Λ is the averagenumber of signals in-band at any given moment (this isreferred to as the duty cycle). However, this was derivedby assuming that the duration of a signal at frequency ν isroughly 1 /ν , which is a good approximation for burstlikesignals, but is inaccurate for the chirp signals emitted bycoalescing compact binaries. As mentioned in Ref. [20],the duration of a compact binary signal in some smallfrequency interval [ ν, ν + δ ν ] is given by∆ t ≈ (cid:18) π / M / ν / (cid:19) − δ ν, (45)where M is the chirp mass of the binary, as defined inEq. (18). Taking this into account, the appropriate limiton the observing time becomes T (cid:29) N pix R + N R ∆ t , (46)where R is the average rate of arrival of GW signals.If we focus on the AGWB at a single frequency, then weare choosing a frequency interval δ ν equal to the frequencyresolution 1 /T . The above implies that the backgroundis only a GRF at this frequency if965 π / M / ν / N R (cid:28) , (47)which is impossible for the sources considered here. It istherefore inconsistent to assume Gaussianity when con-sidering the AGWB at a single frequency. (We note inpassing that due to the large exponent on the frequencyin the equation above, this requirement is much easier tofulfil in the LISA frequency band.)There are two possible ways of addressing this: (i)we can choose a larger frequency interval, integrate thesignal over this interval, and treat the result as a GRF;(ii) we can characterize the anisotropies at a single fre-quency by computing higher-order correlators as well asthe 2PCF (e.g., the AGWB bispectrum and trispectrum).The former is computationally very expensive, due to therequired sampling of many points in GW frequency space(however, it is worth noting that LIGO/Virgo stochasticsearches typically integrate the data over a frequency binthat is significantly larger than 1 /T anyway, so it may de-sirable to compute predictions for this approach regardlessof considerations about Gaussianity). The latter requires a more detailed study, but has the potential to revealmore detailed astrophysical and cosmological information.This will be investigated in a future work.For the purposes of this work, we focus on the 2PCFcomputed at a single frequency, as this still contains agreat deal of important information, and is the first stepin either of the approaches described above. III. ANALYTICAL APPROACH
In this section, we use simple analytical functions for thegalaxy number density ¯ n and galaxy-galaxy 2PCF (cid:104) δ n δ n (cid:105) .This allows us to derive a simple closed-form expressionfor each of the C (cid:96) ’s, up to a frequency-dependent factor A gw which can be integrated numerically. A. Galaxy distribution
In Sec. II we have not yet specified the mean num-ber density distribution ¯ n . For this initial analyticalapproach, we assume that the two galaxy parameters areindependent—i.e., that there is no correlation betweenthe metallicity of a galaxy and its delayed SFR. Thisallows us to write¯ n ( z, ψ d , Z ) = ¯ N ( z ) p ( ψ d | z ) p ( Z| z ) , (48)where p ( ψ d | z ) and p ( Z| z ) are redshift-dependent probab-ility distributions for the two parameters, and¯ N ( z ) ≡ (cid:90) d ζ g ¯ n (49)is the total number density of galaxies per comovingvolume. With this assumption, the density parameter Ω gw reads Ω gw = (cid:88) i π R (local) i I i ( t H ν o ) × (cid:90) z max d z zE ( z ) ¯ N ( z ) (cid:18)(cid:90) ∞ d ψ d ,i p ( ψ d ,i | z ) ψ d ,i (cid:19) × (cid:90) Z max −∞ d Z p ( Z| z ) (cid:90) d ζ b p i ( ζ b ) f Z S i (1 + δ n + ˆ e o · v o ) . (50)Notice that since the integrand is proportional to ψ d ,we only need the first moment (i.e. the mean) of thecorresponding probability distribution,¯ ψ d ≡ (cid:90) ∞ d ψ d p ( ψ d | z ) ψ d , (51)and there is no need to specify anything else about thedistribution. The quantity in Eq. (51) is just the meandelayed SFR, averaged across all galaxies at the appro-priate redshifts. Equivalently, we can think of this as thedelayed form of the cosmic mean SFR per galaxy ¯ ψ ,¯ ψ d = 1ln ( t ( z ) /t min ) (cid:90) t ( z ) t min d(ln t d ) ¯ ψ ( z f ) . (52)Expressions for the cosmic mean SFR are more appro-priately given in units per unit comoving volume ratherthan per galaxy, so we define¯ ψ ( V ) ( z ) ≡ ¯ N ( z ) ¯ ψ ( z ) , (53)with the mean number of galaxies per comoving volume¯ N converting between the two. This function is given inRef. [36] by¯ ψ ( V ) ( z ) = ¯ ψ ( V )peak α exp [ β ( z − z peak )] α − β + β exp [ α ( z − z peak )] , (54)with dimensionless constants α = 2 . β = 2 .
62, andwith the normalization relative to the peak value ¯ ψ ( V )peak =0 . M (cid:12) yr − Mpc − at redshift z peak = 1 .
86. We thuswrite the mean of the delayed SFR distribution as¯ ψ d = 1ln ( t ( z ) /t min ) (cid:90) t ( z ) t min d(ln t d ) ¯ ψ ( V ) ( z f )¯ N ( z f ) . (55)No further information or assumptions about the SFRdistribution are needed. We do, however, need an expres-sion for the total galaxy number density ¯ N ( z ), in orderto convert the SFR per comoving volume ¯ ψ ( V ) back intoSFR per galaxy. We use Ref. [37], where a fit to obser-vational data over a redshift range 0 ≤ z ≤ (cid:18) ¯ N ( z )1 Mpc − (cid:19) = − . − .
08 log (cid:18) t ( z )1 Gyr (cid:19) , (56)where t ( z ) is the age of the Universe at redshift z , as givenby Eq. (24). The dependence of the integrand in Eq. (50) on themetallicity is more complicated due to the factor f Z ; wemust therefore specify the full distribution p ( Z| z ), and notjust its mean. Following Ref. [22], we take the metallicitydistribution as a Gaussian with variance 1 / p ( Z| z ) ∝ exp (cid:104) − (cid:0) Z − ¯ Z ( z ) (cid:1) (cid:105) . (57)Note that this is not strictly a Gaussian, as Z has a finitemaximum value of Z max . The distribution Eq. (57) mustbe reweighted accordingly to give p ( Z| z ) = (cid:113) π exp (cid:104) − (cid:0) Z − ¯ Z ( z ) (cid:1) (cid:105) − erf (cid:2) √ (cid:0) ¯ Z ( z ) − Z max (cid:1)(cid:3) . (58)The cosmic mean metallicity is well modeled by [38]¯ Z ( z ) ≡ log ¯ Z ( z ) Z (cid:12) = Z max + log (cid:32) ψ ( V )norm (cid:90) z max z d z (cid:48) ¯ ψ ( V ) ( z (cid:48) )(1 + z (cid:48) ) E ( z (cid:48) ) (cid:33) , (59)with ¯ ψ ( V ) being the mean SFR from Eq. (54) and ψ ( V )norm anormalizing constant with units of M (cid:12) yr − Mpc − , givenby [38] ψ ( V )norm = ρ bary H o √ y (1 − R ) = 8 . M (cid:12) yr − Mpc − , (60)where ρ bary is the present-day baryon density, y = 0 . R = 0 .
27 is the fraction of metals returned to theinterstellar medium by stars.Given these distributions, the isotropic GW energydensity is therefore¯ Ω gw = (cid:88) i π R (local) i I i ( t H ν o ) (cid:90) z max d z (1 + z ) E ( z ) ln ( t ( z ) /t min ,i ) × (cid:32)(cid:90) t ( z ) t min ,i d(ln t d ) ¯ ψ ( V ) ( z f ) ¯ N ( z )¯ N ( z f ) (cid:33) (cid:90) Z max −∞ d Z p ( Z| z ) (cid:90) d ζ b p i ( ζ b ) f Z S i . (61)with the rate normalizing factor given by I i = (cid:34) t ( z ) /t min ,i ) (cid:32)(cid:90) t ( z ) t min ,i d(ln t d ) ¯ ψ ( V ) ( z f ) ¯ N ( z )¯ N ( z f ) (cid:33) (cid:90) Z max −∞ d Z p ( Z| z ) (cid:90) d ζ b p i ( ζ b ) f Z (cid:35)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z =0 . (62)The result is shown in Fig. 1. ν o / Hz10 − − − − − − π ¯ Ω g w ( ν o ) BNSBBHBHNSTotal (no BHNS)Total (maximum BHNS)
Figure 1. The total GW flux spectrum 4 π ¯ Ω gw given by numerically integrating Eq. (61). (The factor of 4 π is included to givethe total flux, rather than the average flux per solid angle.) The total flux with (without) BHNS mergers is given in black (blue),while the flux coming from BNS, BBH, and BHNS mergers is indicated in red, green, and magenta, respectively. The curvesincluding the BHNS contribution are dashed to indicate that this is an upper bound only. B. Anisotropies
We now consider anisotropies in the AGWB. Since we do not know the form of δ n as a function of sky position,we resort to a statistical description. If the signal is Gaussian, then a complete and simple description is given bythe 2PCF. Even if the signal is non-Gaussian (as suggested by Sec. II C), this is still a convenient and informativedescriptor of the anisotropies. We write the 2PCF of the stochastic background as C gw ( θ o , ν o ) = π ( t H ν o ) Ω (cid:90) z max d z zE ( z ) (cid:90) z max d z (cid:48) z (cid:48) E ( z (cid:48) ) (cid:32)(cid:88) i (cid:90) d ζ g ¯ n (cid:0) z, ζ g (cid:1) (cid:90) d ζ b R i (cid:0) z, ζ g , ζ b (cid:1) S i ( ν s , ζ b ) (cid:33) × (cid:88) j (cid:90) d ζ (cid:48) g ¯ n (cid:0) z (cid:48) , ζ (cid:48) g (cid:1) (cid:90) d ζ (cid:48) b R j (cid:0) z (cid:48) , ζ (cid:48) g , ζ (cid:48) b (cid:1) S j (cid:0) ν (cid:48) s , ζ (cid:48) b (cid:1) (cid:10) δ n (cid:0) z, ˆ e o , ζ g (cid:1) δ n (cid:0) z (cid:48) , ˆ e (cid:48) o , ζ (cid:48) g (cid:1)(cid:11) , (63)which is determined by the final factor on the rhs, namely, the galaxy-galaxy 2PCF, ξ gg ≡ (cid:104) δ n δ n (cid:105) . There exists a0simple analytical estimate for this, ξ gg ( d ) ≈ (cid:18) dd (cid:19) − γ δ ( z − z (cid:48) ) , (64)where d ( z, ˆ e o , ˆ e (cid:48) o ) is the comoving distance between thetwo points, d is the comoving scale at which the cor-relation is unity, and γ is some positive real number.We assume for simplicity that all galaxies cluster inthe same way, regardless of their SFR and metallicity,and that d and γ are constant, and take the values d = (4 . ± . h − Mpc and γ = 1 . ± .
04. Writingthe comoving distance d as d = 2 r tan θ o , θ o ≡ cos − ( ˆ e o · ˆ e (cid:48) o ) , (65)where r = (cid:82) d z /H is the comoving distance between theobserver and the galaxy, the 2PCF can be simply writtenas C gw ( θ o , ν o ) = A gw ( ν o ) (cid:18) tan θ o (cid:19) − γ , (66)where the “amplitude” of the correlation A gw is a functionof the GW frequency alone, and has no angular depend-ence, A gw ( ν o ) ≡ ¯ Ω − (cid:18) d (cid:19) γ (cid:90) z max d zr γ (cid:18) ∂ ¯ Ω gw ∂z (cid:19) . (67)The multipole components C (cid:96) are then given by C (cid:96) ( ν o ) = 2 π A gw (cid:90) +1 − d(cos θ o ) P (cid:96) (cos θ o ) (cid:18) tan θ o (cid:19) − γ , (68)which can be evaluated explicitly (see Appendix A) togive C (cid:96) = 4 π A gw 3 F (cid:0) − (cid:96), (cid:96) + 1 , − γ ; 1 ,
2; 1 (cid:1) sinc( π γ/ , (69)where sinc( x ) ≡ sin ( x ) /x , and F is a generalized hyper-geometric function, which is simple to evaluate numeric-ally.We emphasize that Eq. (69) is only a simple approx-imation of the true angular spectrum of the anisotropies.In reality, the galaxy clustering is more complicated thanwe have assumed, since it evolves with redshift and variesbetween different populations of galaxies with differentphysical attributes. We also note that we have extendedEq. (64) beyond its realm of validity by assuming that itholds for all distances d . The galaxy 2PCF drops below The quoted values are taken from the VIMOS Public ExtragalacticRedshift Survey (VIPERS) [39], and are valid for the brightestgalaxies in a redshift bin centered on z = 1, which are preciselythe galaxies that contribute most strongly to Ω gw . this power law for distances smaller than ≈ . h − Mpcand tends to some finite value as d →
0, and also drops be-low the power law for distances larger than ≈ h − Mpc,eventually becoming negative for large separations. Thereare two resulting inaccuracies in Eqs. (66) and (69) thatare immediately evident. First, the correlation C gw di-verges as θ o →
0, even though this should give a well-defined finite value equal to the variance of the field ateach point, var( δ (s)gw ). Secondly, the (cid:96) = 0 moment ofthe 2PCF evaluates to C = 4 π A gw / sinc( π γ/ >
0, eventhough C should be exactly equal to 0 (as shown inAppendix A). These inaccuracies are due to the extra-polation of Eq. (64) down to d = 0 and out to d → ∞ ,respectively. We therefore turn to a more detailed andaccurate approach in the following section, to address thedeficiencies of this simple model. IV. CATALOGUE APPROACH
The study of anisotropies in the SGWB induced by thelarge-scale distribution of astrophysical sources can bealso performed using a large enough realistic mock cata-logue of galaxies with recipes to infer the production andmerger rates of compact sources. In such catalogues, thegrowth of dark matter structure is first simulated and thehistory of the hierarchical mass assembly is then recordedby means of merger trees of haloes and subhaloes withinsnapshots stored at different time steps. To simulate thevisible galaxies without performing a cosmological hydro-dynamic simulation, one uses a semianalytic model aimingto reproduce the observed properties of these galaxies (e.g.,clustering, counts, scaling relations) at different redshifts.Such models (see, e.g., Refs. [42, 45, 46]) differ by thenumber of astrophysical processes included to describethe fate of baryons within the dark matter haloes, thephenomenological description of these processes, and theirconsistency with the various observational data sets. Inthat respect, quantitative results strongly depend on thetreatment and implementation of feedback mechanisms.Mock light cones can then be generated from the post-processed snapshots of the dynamical collisionless N-bodysimulation, taking care to shuffle them in order to suppressreplication effects.In what follows, we make use of a mock light conecatalogue [41] constructed by applying the “L-galaxies”model [42] to the Millennium simulation [43, 44]. Veryrelevant for our study are the all-sky coverage and depth ofthis mock catalogue, even if the random tiling process usedto generate a mock light cone (the size of which is largerthan the simulated volume) introduces discontinuities inthe density field and hence a loss of clustering information.Hence, the clustering signal suffers from a small negativebias of less than 10% from 1 to 10 h − Mpc around thetwo-point correlation length. In addition, it exhibits finitevolume effects due to the limited size of the simulated boxwhich for instance erases the two-point angular correlationsignal at scales larger than a tenth of the simulation box1 ν o / Hz10 − − − − A g w ( ν o ) BNS onlyBBH onlyBHNS onlyCombined ν o / Hz10 − − − − ¯ Ω w ( ν o ) A g w ( ν o ) Figure 2. The amplitude factor A gw defined in Eq. (67), as a function of GW frequency ν o . The lines in red, green, and magentashow what this factor would be if the background included only BNS, only BBH, or only BHNS events, respectively. The blueline shows the combined result, based on the spectra in Fig. 1 (note that the combined curve is not simply the sum of the BBH,BNS and BHNS curves, due to the normalization with respect to ¯ Ω , which is different in all four cases). The factor diverges atlarge frequencies because it is normalized with respect to the monopole ¯ Ω , which vanishes at these frequencies. The enclosedplot shows ¯ Ω A gw , which encapsulates the absolute size of the anisotropies rather than their relative size compared to themonopole, and therefore does not diverge. size. However, both effects are not important for ourstudy since only a fraction of the spatial information islost affecting scales that are not very relevant for theanisotropies we are looking at. The comoving box ofsize 500 h − Mpc on a side of the Millennium simulationhas indeed been chosen that large in order to encompasstypical scales related to the large-scale features of thespatial distribution of galaxies.The all-sky mock light cone catalogue [41, 42] we usecontains 5,715,694 galaxies, all limited at apparent ABmagnitude of 18 in the r filter from SDSS. It has been builtusing a random tiling technique applied to 64 postpro-cessed snapshots saved during the Millennium simulationwith a time step of about 100 Myr. We queried the data-base and retrieved for each galaxy the following data: its sky location, cosmological redshift, metallicity, andpeculiar velocity. In order to calculate the delayed starformation rate of each galaxy on the light cone, it wasnecessary to access information about its star formationrate at earlier snapshots by querying the full Millenniumsimulation. This was made more complicated by the factthat the light cone galaxies are the result of a sequenceof mergers of smaller galaxies, each with their own in-dependent star formation history. In order to accountfor this, we queried the Millennium simulation to extractthe full star formation history of each light cone galaxy,given by the star formation rates and redshifts of its pro-genitors at earlier snapshots. This included a total of973,224,532 redshift and SFR measurements from theprogenitor galaxies.2
GW energy overdensity δ (s)gw -1.00e+00 1.00e+00 Figure 3. A HEALPix [40] map of the GW energy overdensity δ (s)gw constructed from the all-sky mock light cone catalogue [41–44],as described in Sec. IV. This was generated with the HEALPix N side parameter set to 256, corresponding to an angular resolutionof 13.7 arcminutes, and an average of 7.3 galaxies per pixel. We note that due to the magnitude-limited sample usedto construct the light cone catalogue, it only extends outto a redshift of z ≈ .
78. While this is sufficient to studyanisotropies on the angular scales we are interested in,it means that the total energy density ¯ Ω gw given by thecatalogue will be significantly less than the true value,due to the missing contribution from redshifts z > . δ (s)gw , this will have minimal effecton the C (cid:96) spectrum. The main advantage of the mockcatalogue lies in its accurate nonlinear modeling of thegalaxy clustering statistics, which remains valid on thescales of interest even when the redshift limit is introduced.In performing our analysis of the mock catalogue, weadhere to the (now outdated) WMAP values of the cosmo-logical parameters that were used in the Millennium sim-ulation: H o = 73 km s − Mpc − , Ω m = 0 . Ω Λ = 0 . γ and d that match those of the simula-tion, which are themselves consistent with the 2-degreeField Galaxy Redshift Survey [43, 47]: γ = 1 . ± . d = 5 . ± . h − Mpc. This ensures consistencybetween our results and the underlying simulation. Weconfirm that replacing these with the more up-to-date values mentioned previously ( H o = 67 . − Mpc − , Ω m = 0 . Ω Λ = 0 . d = 4 . h − Mpc, γ = 1 . A. Reconstructing the SGWB from pointlikesources
In order to use the information extracted from thecatalogue, we must first explicitly rewrite our equationsfor Ω gw in terms of the available data for each galaxy. Wedo this by expressing the galaxy number density n as aweighted sum of Dirac delta functions, reflecting the factthat we treat each galaxy as a point source.Let us consider a catalogue of N galaxies, indexed by alabel k . We write their redshift, sky location, delayed starformation rate, and log-normalized metallicity as z k , ˆ e k ,¯ ψ d ,k ( z ), and Z k , respectively. We also allow each galaxyto have a peculiar velocity v k along the line of sight. Thismeans that galaxy k has a source-frame frequency givenby ν s ,k = ν o (1 + z k )(1 + v k − ˆ e k · v o ) . (70)By integrating the number density per comoving volume3 n over redshift, SFR, and metallicity, we must have N = (cid:90) z max z =0 d V ( z ) (cid:90) d ζ g n = (cid:90) S d σ o (cid:90) z max d z r H (cid:90) ∞ d ¯ ψ d (cid:90) Z max −∞ d Z n. (71) This implies that n = H o (cid:88) k E ( z k ) r ( z k ) δ (2) ( ˆ e o , ˆ e k ) δ ( z − z k ) δ (cid:0) ¯ ψ d − ¯ ψ d ,k ( z k ) (cid:1) δ ( Z − Z k ) . (72)The SGWB is then given by Ω gw = (cid:88) i π t H ν o ) (cid:90) z max d z zE ( z ) (cid:90) d ζ g n (1 + ˆ e o · v o ) (cid:90) d ζ b R i ( z, Z , ζ b ) S i ( ν s , ζ b )= (cid:88) k (cid:88) i π H o t H ν o ) z k r ( z k ) (1 + ˆ e k · v o ) (cid:90) d ζ b R i ( z k , Z k , ζ b ) S i ( ν s ,k , ζ b ) δ (2) ( ˆ e o , ˆ e k ) . (73)Let us remind the reader that i ∈ { BNS , BBH , BHNS } indexes a sum over the different types of binary merger ineach galaxy, while k ∈ { , , . . . , N } indexes a sum over the galaxies in the catalogue. The above expression for Ω gw becomes much simpler if we define the weight w k ( ν o ) ≡ (cid:88) i π H o t H ν o ) z k r ( z k ) (1 + ˆ e k · v o ) (cid:90) d ζ b R i ( z k , Z k , ζ b ) S i ( ν s ,k , ζ b ) (74)for each galaxy, so that Ω gw = (cid:88) k w k δ (2) ( ˆ e o , ˆ e k ) . (75)In the analytic approach we took in the previous section,we were limited by the fact that δ n was not known as afunction of sky location; all we knew was the (approxim-ate) 2PCF of δ n , so we were forced to use the 2PCF of Ω gw . Using a catalogue we now have an explicit expres-sion for Ω gw as a function of sky location, and we can usewhatever statistics we like to describe it.One simple choice is to calculate the spherical multipolecomponents of Ω gw , Ω (cid:96)m ( ν o ) ≡ (cid:90) S d σ o Ω gw Y ∗ (cid:96)m ( ˆ e o ) = (cid:88) k w k Y ∗ (cid:96)m ( ˆ e k ) . (76)However, these quantities become expensive to computeif our catalogue contains a large number of galaxies, aseach component requires N evaluations of a sphericalharmonic function. This work focuses on the 2PCF of theGW energy overdensity C gw = (cid:68) δ (s)gw δ (s)gw (cid:69) , which capturesmuch of the information contained in the anisotropies.It is possible to write down an expression for the C (cid:96) components of this function in terms of the weights w k and position vectors ˆ e k of each galaxy; however, these areeven more expensive to compute, as they require a sumover all pairs of galaxies which scales as N − N , or morethan 3 × operations for each C (cid:96) . In practice, it is much more efficient to compute the stat-istics of the background using HEALPix [40], a powerfulsoftware package for manipulating pixelized maps of thesphere. We construct a HEALPix map by computing theweight w k of each galaxy [given by Eq. (74)] and addingthis to the pixel corresponding to the galaxy’s sky position ˆ e k . The C (cid:96) ’s (or other relevant quantities) can then becalculated directly in a matter of seconds using HEALPixroutines. V. RESULTS AND DISCUSSION
We have numerically integrated Eq. (61) to find thefrequency spectrum of the isotropic component of theAGWB as shown in Fig. 1, which matches the correspond-ing figure in Ref. [22]. At frequencies up to ≈
20 Hz, allsources are included and are in the inspiral regime, givinga monopole ¯ Ω gw that scales as ν / . The number of emit-ting sources decreases at higher frequencies, causing thespectrum to taper off and eventually vanish at ≈ A gw given in Eq. (67), and hence obtain a simple ana-lytical prediction for the C (cid:96) spectrum of the backgroundusing Eq. (69). The value of this factor as a function of http://healpix.sourceforge.net ‘ − − − − ‘ ( ‘ + ) C ‘ / π CatalogueAnalytical
Figure 4. The quantity (cid:96) ( (cid:96) + 1) C (cid:96) / π , which is the approximate contribution to the anisotropic variance of δ (s)gw as a function ofln (cid:96) . The red curve shows the simple analytical prediction given in Eq. (69), while the blue curve shows the spectrum computedfrom the map in Fig. 3 using HEALPix. Both curves include error regions from cosmic variance, while the blue curve includesPoisson errors associated with the finite number of galaxies per pixel in the HEALPix map. frequency is shown in Fig. 2. In the inspiral-dominatedregime below ≈
20 Hz, the value of A gw is essentially con-stant at ≈ × − . As the frequency increases above thisregime, the falloff in the number of sources increases thegranularity of the background, and therefore the size ofthe anisotropies, causing A gw to increase by many ordersof magnitude, and eventually diverge as ¯ Ω gw → C (cid:96) at any fre-quency in the LIGO-Virgo band. For concreteness, wefocus on a frequency of 65 .
75 Hz, as this has been used asthe reference frequency in previous directional searchesfor an inspiral-dominated astrophysical background con-ducted by LIGO/Virgo [17]. We find A gw ( ν o = 65 .
75 Hz) ≈ . × − . (77)This includes BHNS mergers at the maximal rate given inEq. (29); setting the BHNS rate to 0 gives a 6% decreasein A gw . We have also computed the C (cid:96) spectrum of the AGWBfrom an all-sky mock light cone catalogue based on theMillennium simulation, giving a much more accurate de-scription of the anisotropies by relaxing many of the as-sumptions about the galaxy-galaxy 2PCF that were usedto derive Eq. (61). This was done by constructing a GWoverdensity map as described in Sec. IV, which is shownin Fig. 3. The C (cid:96) spectrum computed from this map usingthe HEALPix software package is shown in Fig. 4, alongwith the analytical prediction. Despite the simplicity ofthe analytical model, and the fact that it contains nofree parameters, it is in very strong agreement with thecorresponding result from the catalogue at the largestangular scales (i.e., the lowest (cid:96) -modes). It also capturesthe scaling of C (cid:96) with (cid:96) at small scales (i.e., large (cid:96) ), al-though it overestimates the anisotropic variance at thesescales by a constant factor of ≈
4. It is worth mentioningthat such high values of (cid:96) are inaccessible to the current5detector network due to their poor angular resolution—current directional searches only probe the first few (cid:96) [17],although this is expected to improve as more detectorsare added to the network. Note that in order to ensure afair comparison in Fig. 4, the analytical curve shown isfor the cosmological parameters H o , Ω m , Ω Λ and galaxy2PCF parameters γ , d corresponding to the simulation.Replacing these with the more up-to-date values fromPlanck [48] and VIPERS [39] has negligible effect.The results in Fig. 4 can be compared directly withFig. 2 of Ref. [21] [with the caveat that Fig. 4 includesBNS and BHNS mergers as well as BBH, whereas Ref. [21]includes only BBH—however, Fig. 2 indicates that thisshould only introduce a factor of O (1) between the two].We note that the C (cid:96) ’s in Fig. 4 are larger in amplitudethan in Ref. [21], and do not fall off as quickly with (cid:96) .This perhaps indicates that the linear perturbation-theoryapproach adopted in Ref. [21] to describe the galaxy dis-tribution underestimates the clustering at small scales,where the perturbations to the density field become nonlin-ear. A more detailed comparison of the two sets of resultsis needed to fully assess the accuracy and applicability ofeach approach.Figure 4 can also be compared with Fig. 5 of Ref. [20],which shows the corresponding C (cid:96) spectra for cosmicstring networks with different string tensions Gµ . Theanisotropies in the astrophysical background are clearlymuch larger in amplitude, and become nonlinear for mod-erate values of (cid:96) . This is unsurprising, given that the bulkof the GW flux from astrophysical sources is emitted atmuch lower redshifts, so that angular separations on thesky correspond to much smaller physical separations thanin the cosmic string case, and therefore to scales wherethe departure from homogeneity is more evident.We have also computed the kinematic dipole factor,which in the case shown in Figs. 3 and 4 has a value of D ( ν o = 65 .
75 Hz) ≈ . × − . (78)This is an order of magnitude smaller than in the cosmicstring case, which is intuitively sensible given that thegalaxies we consider are generally at much lower redshiftsand therefore their velocities due to the Hubble flow aremuch smaller, giving a weaker Doppler effect. The factthat D is smaller in this case, while δ (s)gw is simultaneouslylarger, means that the effects of the kinematic dipole areunimportant for the astrophysical background. VI. CONCLUSION
We have developed a detailed anisotropic model forthe AGWB, including the most important sources for theLIGO-Virgo frequency band (BBH, BNS, and BHNS).The angular spectrum of the anisotropies, quantified bythe C (cid:96) components, has been calculated through two com-plementary approaches: a simple, closed-form analyticalexpression Eq. (69), and a detailed numerical study using an all-sky mock light cone galaxy catalogue from the Mil-lennium simulation [41–44]. The two approaches are inexcellent agreement at large angular scales ( (cid:96) (cid:46) (cid:96) ( (cid:96) + 1) C (cid:96) / π over many orders ofmagnitude despite the simplicity of the analytical modeland the lack of free parameters. These anisotropies areconsiderably larger in amplitude than those in the temper-ature of the CMB, or those in the SGWB due to cosmicstrings [20], and become nonlinear at higher multipoles (cid:96) .This shows that modeling the AGWB in a purely isotropicmanner neglects a great deal of astrophysical and cosmo-logical information, thereby motivating future theoreticaland observational work. We expect that in the near fu-ture, anisotropies in the AGWB (and the SGWB moregenerally) will become an important probe of the large-scale structure of the Universe, and of the astrophysicalprocesses that occur within it.We have highlighted several key avenues for future workto explore. One of the key differences between the SGWBand the CMB is the very broad frequency spectrum overwhich the SGWB can be investigated. This motivatesthe development of models that are valid at frequenciesoutside of the LIGO-Virgo band; the study of AGWBanisotropies in the LISA band will be of particular interestin the coming decades. As shown in Sec. II C, the use ofnarrow frequency bins makes it impossible to guaranteethat the AGWB is Gaussian, and a full characterization ofthe anisotropies in this bin requires higher-order correlat-ors of the field (such as the bispectrum and trispectrum).This is a daunting task, both theoretically and observa-tionally, but it holds the promise of incredibly rich newastrophysical and cosmological information.There will also be a need to develop increasingly ac-curate models of the AGWB anisotropies. This will beaddressed in future work by improving upon the fiducialastrophysical model of Ref. [22], using past and futureLIGO-Virgo observing runs to enhance our knowledge ofthe relevant sources. The catalogue approach described inSec. IV will also be improved upon by using more accurateand complete galaxy catalogues. It will also be interestingto explore different ways in which the statistical proper-ties of the AGWB might differ from those of the modeldescribed here; e.g., the extension to models which arestatistically nonstationary, or which are no longer stat-istically isotropic due to the inclusion of sources in thegalactic plane of the Milky Way.Finally, we plan to produce realistic mock data from thismodel and from future models, and use them to test andoptimize data analysis and parameter estimation methods,leading the way to future detections and measurement ofAGWB anisotropies. ACKNOWLEDGMENTS
We thank Joe Romano and Andrew Matas for readingthe manuscript carefully and providing us with valuable6comments. The Millennium Simulation databases usedin this paper and the web application providing onlineaccess to them were constructed as part of the activitiesof the German Astrophysical Virtual Observatory. Someof the results in this paper have been derived using theHEALPix package [40]. A.C.J. is supported by King’sCollege London through a Graduate Teaching Scholarship.M.S. is supported in part by the Science and TechnologyFacility Council (STFC), United Kingdom, under GrantNo. ST/P000258/1.
Appendix A: Explicitly evaluating the analytical C (cid:96) spectrum The C (cid:96) ’s are given by C (cid:96) ( ν o ) = 2 π A gw (cid:90) +1 − d(cos θ o ) P (cid:96) (cos θ o ) (cid:18) tan θ o (cid:19) − γ . (A1)Applying simple trigonometry, we havetan θ o (cid:114) − cos θ o θ o , and therefore C (cid:96) = 2 π A gw (cid:90) +1 − d x P (cid:96) ( x ) (cid:18) x − x (cid:19) γ/ . This can be evaluated by writing the Legendre polynomialas a sum, P (cid:96) ( x ) = (cid:96) (cid:88) k =0 ( (cid:96) + k )!( k !) ( (cid:96) − k )! (cid:18) x − (cid:19) k , so that we have C (cid:96) = 2 π A gw (cid:96) (cid:88) k =0 ( (cid:96) + k )!( k !) ( (cid:96) − k )! (cid:90) +1 − d x (cid:18) x − x (cid:19) γ/ (cid:18) x − (cid:19) k . With a change of variables to y ≡ − x , the integral be-comes (cid:90) +1 − d x (cid:18) x − x (cid:19) γ/ (cid:18) x − (cid:19) k = 2( − k (cid:90) +10 d y (1 − y ) γ/ y k − γ . This is a Euler integral of the first kind, and evaluates to (cid:90) +10 d y (1 − y ) γ/ y k − γ = Γ (cid:0) k − γ (cid:1) Γ (cid:0) γ (cid:1) Γ ( k + 2) , where we assume 0 < γ <
2. We therefore have C (cid:96) = 4 π A gw (cid:96) (cid:88) k =0 ( − k ( (cid:96) + k )!( k !) ( (cid:96) − k )! Γ (cid:0) k − γ (cid:1) Γ (cid:0) γ (cid:1) Γ ( k + 2)= 4 π A gw (cid:96) (cid:88) k =0 ( − k Γ ( (cid:96) + k + 1) Γ ( k + 1) Γ ( (cid:96) − k + 1) Γ (cid:0) k − γ (cid:1) Γ (cid:0) γ (cid:1) Γ ( k + 2)= 4 π A gw ∞ (cid:88) k =0 ( − k Γ ( (cid:96) + k + 1) Γ ( k + 1) Γ ( (cid:96) − k + 1) Γ (cid:0) k − γ (cid:1) Γ (cid:0) γ (cid:1) Γ ( k + 2) , where the final line follows because Γ ( (cid:96) + k + 1) /Γ ( (cid:96) − k + 1) vanishes for k > (cid:96) . Wecan simplify further by defining the rising Pochammersymbol, ( x ) k ≡ x ( x + 1) · · · ( x + k − , and by using Euler’s reflection formula,sinc( π x ) ≡ sin π x π x = 1 Γ (1 + x ) Γ (1 − x ) . This gives C (cid:96) = 4 π A gw sinc( π γ/ ∞ (cid:88) k =0 ( − (cid:96) ) k ( (cid:96) + 1) k (cid:0) − γ (cid:1) k (1) k (2) k k ! . The series above defines a generalized hypergeometricfunction F , so that we get our final expression, C (cid:96) = 4 π A gw 3 F (cid:0) − (cid:96), (cid:96) + 1 , − γ ; 1 ,
2; 1 (cid:1) sinc( π γ/ . (A2)Note that when (cid:96) = 0, the hypergeometric functionevaluates to unity, and we have C = 4 π A gw sinc( π γ/ > . (A3)However, the C component is just an average of C gw over7the sphere, C = 2 π (cid:90) +1 − d(cos θ o ) C gw = (cid:90) π d φ o (cid:90) π d θ o sin θ o C gw = (cid:90) S d σ o C gw , and therefore must be equal to 0, since the sphericalintegral and the ensemble averaging process commute, andsince the average of δ (s)gw over the sphere is 0 by definition, C = (cid:28) δ (s)gw (cid:90) S d σ o δ (s)gw (cid:29) = 0 . (A4) This is indicative of the inaccuracies inherent to Eq. (69). Appendix B: Waveform numerical constants
The following expressions and numerical values areneeded to fully specify the hybrid waveforms given inEq. (17). They exactly match those given in Ref. [32]. i ν i y (10) i y (11) i y (12) i y (20) i y (21) i y (30) i − . − χ ) . + 3 . − χ ) . +0 . . − . − . − . − . − . − χ ) . ] / . − . − . − . . . − . − χ ) . ](1 − χ ) . / − . − . . . − . − .
874 0 . . χ + 0 . χ − . − . . − . . . α = − (cid:18) M M (cid:19) / ,(cid:15) = 1 . χ − . ,c = ν − (cid:34) (cid:80) i =2 α i ( π GM ν ) i/ (cid:80) i =1 (cid:15) i ( π GM ν ) i/ (cid:35) , α = (cid:34) − (cid:18) M M (cid:19) / (cid:35) χ,(cid:15) = − . χ + 1 . c = c ν − / (cid:34) (cid:88) i =1 (cid:15) i ( π GM ν ) i/ (cid:35) . [1] G. M. Harry et al. (LIGO Scientific), Class. Quant. Grav. , 084006 (2010).[2] J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. ,074001 (2015), arXiv:1411.4547 [gr-qc].[3] F. Acernese et al. (VIRGO), Class. Quant. Grav. ,024001 (2015), arXiv:1408.3978 [gr-qc].[4] B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev.Lett. , 241103 (2016), arXiv:1606.04855 [gr-qc].[5] B. P. Abbott et al. (VIRGO, LIGO Scientific), Phys. Rev.Lett. , 221101 (2017), arXiv:1706.01812 [gr-qc].[6] B. P. Abbott et al. (Virgo, LIGO Scientific), Astrophys.J. , L35 (2017), arXiv:1711.05578 [astro-ph.HE].[7] B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev.Lett. , 141101 (2017), arXiv:1709.09660 [gr-qc].[8] B. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev.Lett. , 161101 (2017), arXiv:1710.05832 [gr-qc].[9] S. W. Ballmer, Class. Quant. Grav. , S179 (2006),arXiv:gr-qc/0510096 [gr-qc].[10] B. Abbott et al. (LIGO Scientific), Phys. Rev. D76 ,082003 (2007), arXiv:astro-ph/0703234 [ASTRO-PH].[11] E. Thrane, S. Ballmer, J. D. Romano, S. Mitra, D. Taluk-der, S. Bose, and V. Mandic, Phys. Rev.
D80 , 122002(2009), arXiv:0910.0858 [astro-ph.IM].[12] J. Abadie et al. (LIGO Scientific), Phys. Rev. Lett. ,271102 (2011), arXiv:1109.1809 [astro-ph.CO].[13] C. Messenger, H. J. Bulten, S. G. Crowder, V. Dergachev, D. K. Galloway, E. Goetz, R. J. G. Jonker, P. D. Lasky,G. D. Meadors, A. Melatos, S. Premachandra, K. Riles,L. Sammut, E. H. Thrane, J. T. Whelan, and Y. Zhang,Phys. Rev.
D92 , 023006 (2015), arXiv:1504.05889 [gr-qc].[14] C. Chung, A. Melatos, B. Krishnan, and J. T. Whelan,Mon. Not. Roy. Astron. Soc. , 2650 (2011),arXiv:1102.4654 [gr-qc].[15] L. Sun, A. Melatos, P. D. Lasky, C. T. Y. Chung,and N. S. Darman, Phys. Rev.
D94 , 082004 (2016),arXiv:1610.00059 [gr-qc].[16] J. Aasi et al. (VIRGO, LIGO Scientific), Phys. Rev.
D88 ,122004 (2013), arXiv:1309.6160 [astro-ph.HE].[17] B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev.Lett. , 121102 (2017), arXiv:1612.02030 [gr-qc].[18] G. Cusin, C. Pitrou, and J.-P. Uzan, Phys. Rev.
D96 ,103019 (2017), arXiv:1704.06184 [astro-ph.CO].[19] G. Cusin, C. Pitrou, and J.-P. Uzan, Phys. Rev.
D97 ,123527 (2018), arXiv:1711.11345 [astro-ph.CO].[20] A. C. Jenkins and M. Sakellariadou, (2018),arXiv:1802.06046 [astro-ph.CO].[21] G. Cusin, I. Dvorkin, C. Pitrou, and J.-P. Uzan, Phys.Rev. Lett. , 231101 (2018), arXiv:1803.03236 [astro-ph.CO].[22] B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev.Lett. , 091101 (2018), arXiv:1710.05837 [gr-qc]. [23] I. Dvorkin, J.-P. Uzan, E. Vangioni, and J. Silk, (2017),10.1093/mnras/sty1414, arXiv:1709.09197 [astro-ph.HE].[24] I. Dvorkin, J.-P. Uzan, E. Vangioni, and J. Silk,Phys. Rev. D94 , 103011 (2016), arXiv:1607.06818 [astro-ph.HE].[25] I. Dvorkin, E. Vangioni, J. Silk, J.-P. Uzan, and K. A.Olive, Mon. Not. Roy. Astron. Soc. , 3877 (2016),arXiv:1604.04288 [astro-ph.HE].[26] P. Amaro-Seoane et al. (LISA), (2017), arXiv:1702.00786[astro-ph.IM].[27] T. Regimbau, M. Evans, N. Christensen, E. Katsavounidis,B. Sathyaprakash, and S. Vitale, Phys. Rev. Lett. ,151105 (2017), arXiv:1611.08943 [astro-ph.CO].[28] R. Smith and E. Thrane, Phys. Rev. X8 , 021019 (2018),arXiv:1712.00688 [gr-qc].[29] T. Namikawa, A. Nishizawa, and A. Taruya, Phys.Rev. Lett. , 121302 (2016), arXiv:1511.04638 [astro-ph.CO].[30] T. Namikawa, A. Nishizawa, and A. Taruya, Phys. Rev. D94 , 024013 (2016), arXiv:1603.08072 [astro-ph.CO].[31] P. Ajith, S. Babak, Y. Chen, M. Hewitson, B. Krishnan,A. M. Sintes, J. T. Whelan, B. Brugmann, P. Diener,N. Dorband, J. Gonzalez, M. Hannam, S. Husa, D. Pollney,L. Rezzolla, L. Santamaria, U. Sperhake, and J. Thorn-burg, Phys. Rev.
D77 , 104017 (2008), [Erratum: Phys.Rev.D79,129901(2009)], arXiv:0710.2335 [gr-qc].[32] P. Ajith, M. Hannam, S. Husa, Y. Chen, B. Brugmann,N. Dorband, D. Muller, F. Ohme, D. Pollney, C. Reisswig,L. Santamaria, and J. Seiler, Phys. Rev. Lett. ,241101 (2011), arXiv:0909.2867 [gr-qc].[33] B. P. Abbott et al. (Virgo, LIGO Scientific), (2018),arXiv:1805.11579 [gr-qc].[34] B. P. Abbott et al. (Virgo, LIGO Scientific), Astrophys.J. , L21 (2016), arXiv:1607.07456 [astro-ph.HE].[35] J. Abadie et al. (VIRGO, LIGO Scientific), Class. Quant. Grav. , 173001 (2010), arXiv:1003.2480 [astro-ph.HE].[36] E. Vangioni, K. A. Olive, T. Prestegard, J. Silk,P. Petitjean, and V. Mandic, Mon. Not. Roy. Astron.Soc. , 2575 (2015), arXiv:1409.2462 [astro-ph.GA].[37] C. J. Conselice, A. Wilkinson, K. Duncan, andA. Mortlock, Ap. J. , 83 (2016), arXiv:1607.03909.[38] K. Belczynski, D. E. Holz, T. Bulik, andR. O’Shaughnessy, Nature , 512 (2016),arXiv:1602.04531 [astro-ph.HE].[39] F. Marulli et al. , Astron. Astrophys. , A17 (2013),arXiv:1303.2633 [astro-ph.CO].[40] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wan-delt, F. K. Hansen, M. Reinecke, and M. Bartelmann,Astrophys. J. , 759 (2005), arXiv:astro-ph/0409513[astro-ph].[41] J. Blaizot, Y. Wadadekar, B. Guiderdoni, S. T. Colombi,E. Bertin, F. R. Bouchet, J. E. G. Devriendt, andS. Hatton, Mon. Not. Roy. Astron. Soc. , 159(2005), arXiv:astro-ph/0309305 [astro-ph].[42] G. De Lucia and J. Blaizot, Mon. Not. Roy. Astron. Soc. , 2 (2007), arXiv:astro-ph/0606519 [astro-ph].[43] V. Springel et al. , Nature , 629 (2005), arXiv:astro-ph/0504097 [astro-ph].[44] G. Lemson (VIRGO), (2006), arXiv:astro-ph/0608019[astro-ph].[45] S. Hatton, J. E. G. Devriendt, S. Ninin, F. R. Bouchet,B. Guiderdoni, and D. Vibert, Mon. Not. Roy. Astron.Soc. , 75 (2003), arXiv:astro-ph/0309186 [astro-ph].[46] B. M. B. Henriques, S. White, P. Thomas, R. Angulo,Q. Guo, G. Lemson, V. Springel, and R. Overzier,Mon. Not. Roy. Astron. Soc. , 2663 (2015),arXiv:1410.0365 [astro-ph.GA].[47] E. Hawkins et al. , Mon. Not. Roy. Astron. Soc. , 78(2003), arXiv:astro-ph/0212375 [astro-ph].[48] P. A. R. Ade et al. (Planck), Astron. Astrophys.594