Exploring galaxies-gravitational waves cross-correlations as an astrophysical probe
PPrepared for submission to JCAP
Exploring galaxies-gravitationalwaves cross-correlations as anastrophysical probe
Giulio Scelfo, a,b,c
Lumen Boco, a,b,c
Andrea Lapi, a,b,c,d
MatteoViel a,b,c,d a SISSA, Via Bonomea 265, 34136 Trieste, Italy b INFN, Sezione di Trieste, Via Bonomea 265, 34136 Trieste, Italy c IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2, 34151, Trieste, Italy d INAF/OATS, Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34143 Trieste, ItalyE-mail: [email protected], [email protected], [email protected], [email protected]
Abstract.
Gravitational waves astronomy has opened a new opportunity to study the Uni-verse. Full exploitation of this window can especially be provided by combining data comingfrom gravitational waves experiments with luminous tracers of the Large Scale Structure, likegalaxies. In this work we investigate the cross-correlation signal between gravitational wavesresolved events, as detected by the Einstein Telescope, and actively star-forming galaxies.The galaxies distribution is computed through their UV and IR luminosity functions andthe gravitational waves events, assumed to be of stellar origin, are self-consistently computedfrom the aforementioned galaxies distribution. We provide a state-of-the-art treatment bothon the astrophysical side, taking into account the impact of the star formation and chemicalevolution histories of galaxies, and in computing the cross-correlation signal, for which weinclude lensing and relativistic effects. We find that the measured cross-correlation signalcan be sufficiently strong to overcome the noise and provide a clear signal. As a possibleapplication of this methodology, we consider a proof-of-concept case in which we aim atdiscriminating a metallicity dependence on the compact objects merger efficiency against areference case with no metallicity dependence. When considering galaxies with a Star For-mation Rate ψ > M (cid:12) / yr , a Signal-to-Noise ratio around a value of 2-4 is gained after adecade of observation time, depending on the observed fraction of the sky. This formalismcan be exploited as an astrophysical probe and could potentially allow to test and comparedifferent astrophysical scenarios. a r X i v : . [ a s t r o - ph . C O ] S e p ontents After the first detection of a Gravitational Wave (GW) signal was announced, originatingfrom the merger of Binary Black Holes (BBH) of a total mass M tot ∼ M (cid:12) [1, 2], the era ofGW astronomy began. Its groundbreaking importance comes from the fact that it opened acompletely new way to observe the cosmos, using an observable that could not be exploitedbefore. Among the newly opened directions, the birth of GW astronomy also led to newpossibilities in the multimessenger field of tracers cross-correlations.The study of cross-correlations between distinct tracers is not new. Indeed, severalstudies regarding e.g., correlations between the Large Scale Structure (LSS) and the CosmicMicrowave Background (see e.g., refs. [3–8]), neutrinos (see e.g., ref. [9]) or among differentLSS tracers (see e.g., refs. [10–13]) have been performed. For what concerns the specific pos-sibility to cross-correlate GW signals with LSS tracers, various works on different applicationshave been made, such as the investigation of the origin of merging BBHs [14, 15], the study ofanisotropies of the number density and luminosity distances of compact binaries [16] and ofthe stochastic GW background [17, 18], the investigation of the GW bias [19], the possibilityof alternatives to General Relativity [20] and several others (see e.g., refs. [21, 22]).– 1 –n this paper we extend the work addressed by the community in the GW × LSS areastudying the measurable cross-correlation signal with a refined description of both the LSSand GW tracers. Regarding the characterization of the LSS tracers we make use of activelystar-forming galaxies. Observations of the last decade, with the advent of high redshift far-IR/sub-mm surveys, have helped in robustly characterizing the galaxies luminosity functions,from which the Star Formation Rate (SFR) ψ can be derived, allowing us to have a rather solidstatistics of actively star-forming galaxies (for a more detailed discussion see section 3.1.1 andreferences therein). Therefore, we exploit the galaxy SFR to organize different galaxy typesin different SFR bins. In this way we can look at the contribution to the cross-correlationsignal coming from galaxies with different star formation activities.On the other hand, we model the redshift distributions of the detected GW signals,coming from all types of merging compact objects COs (BH-BH, BH-NS, NS-NS, where NSstands for Neutron Stars), by convolving a detector sensitivity curve with the intrinsic mergingrates. We self-consistently derive the latters from the aforementioned galaxies distribution.Thus, the two tracers considered here are not coming from different and independent sources:we are looking at the same objects (galaxies) but through different messengers (light andGWs). However, we need to take into account that the GWs distribution does not only dependon the SFR of the galaxies. Both the number of merging compact objects and their chirpmass (affecting the GW signal detectability) can strongly depend also on the environmentalconditions in which the binary forms and evolves. For these reasons, in this work we use arefined determination of the COs merging rates, following reference [23], where a metallicitydistribution is associated to each single galaxy through a chemical evolution model. Once theintrinsic merging CO distribution is computed, we convolve it with the sensitivity curve ofthe future third generation GW observatory Einstein Telescope (ET) [24] to get the detectedGW events redshift distribution.These characterizations altogether lead to our forecast on the cross-correlation signalthat can be obtained by realistically modeling these two types of tracers, especially awaitingthe soon-to-come GWs detections from third generation observatories. We perform here botha tomographic and a non tomographic approach.As a possible application of this cross-correlation formalism, we explore the idea of testingdifferent astrophysical scenarios, which predict different GWs distributions, clustering andother specifics. We consider an exemplificative proof-of-concept case, in which we test whethermetallicity dependencies on the COs merging efficiency can be detected and distinguishedwith respect to a benchmark case with no metallicty dependence. Looking at the distributionand clustering of these two tracers can be a promising tool to discriminate between differentfeatures imprinted by several astrophysical mechanisms.This work is structured as follows: in section 2 we provide the mathematical backgroundto describe the cross-correlation signal, given by the number counts angular power spectra C (cid:96) ’s; in section 3 we describe the tracers considered and the theoretical background behindthe physical quantities characterizing them; in section 4 we present estimates on the cross-correlation signals through Signal-to-Noise computations; in section 5 we describe how topotentially distinguish astrophysical features from different models and predict its viabilityusing a test scenario; in section 6 we draw our conclusions.– 2 – Cross-correlation formalism
We describe the cross-correlation between two tracers by working in the harmonic space,considering as observable the number counts angular power spectrum C (cid:96) . The multipolenumber (cid:96) relates to the angular resolution θ as (cid:96) ∼ o /θ (see e.g., refs. [25–35] for worksabout giving up the flat sky approximation and the advantages of working in the harmonicspace). What follows is valid both in the case when tomographic maps of the two tracers areavailable and in the case in which tomography is not performed, whereas the latter case cansimply be seen as the first one reduced to a single redshift bin.Being X and Y the two tracers, we can write the relation between the observed angularpower spectrum ˜ C XY(cid:96) ( z i , z j ) (obtained cross-correlating tracer X in redshift bin z i with tracer Y in bin z j ) and the harmonic coefficients a (cid:96)m as: (cid:104) a X(cid:96)m ( z i ) a Y ∗ (cid:96) (cid:48) m (cid:48) ( z j ) (cid:105) = δ (cid:96)(cid:96) (cid:48) δ mm (cid:48) ˜ C XY(cid:96) ( z i , z j ) , (2.1)where δ is the Kronecker delta. The observed harmonic coefficients a X(cid:96)m ( z i ) are given by thesum of the partial wave coefficients of the signal and of the noise: a X(cid:96)m ( z i ) = s X(cid:96)m ( z i ) + n X(cid:96)m ( z i ) . (2.2)The observed angular power spectra read as ˜ C XY(cid:96) ( z i , z j ) = C XY(cid:96) ( z i , z j ) + δ XY δ ij N X(cid:96) ( z i ) . (2.3)The angular power spectrum is directly obtained from the signal wave coefficients as [36, 37] (cid:104) s X(cid:96)m ( z i ) s Y ∗ (cid:96) (cid:48) m (cid:48) ( z j ) (cid:105) = δ (cid:96)(cid:96) (cid:48) δ mm (cid:48) C XY(cid:96) ( z i , z j ) . (2.4)We construct the noise angular power spectrum from the shot noise N X(cid:96) ( z i ) (inversely pro-portional to the number of sources per steradian), assuming no other sources of error and nocorrelation between noise terms of different experiments and z bins. The expectation valueof the noise can then be written as (cid:104) n X(cid:96)m ( z i ) n Y ∗ (cid:96) (cid:48) m (cid:48) ( z j ) (cid:105) = δ (cid:96)(cid:96) (cid:48) δ mm (cid:48) δ XY δ ij N X(cid:96) ( z i ) . (2.5)Assuming signal and noise as statistically independent, one can write (cid:104) s X(cid:96)m ( z i ) n Y ∗ (cid:96) (cid:48) m (cid:48) ( z j ) (cid:105) = 0 . (2.6)Finally, the angular power spectrum C XY(cid:96) ( z i , z j ) for different tracers and different redshiftbins can be written as C XY(cid:96) ( z i , z j ) = 2 π (cid:90) dkk P ( k )∆ X,z i (cid:96) ( k )∆ Y,z j (cid:96) ( k ) , (2.7)where P ( k ) = k P ( k ) is the primordial power spectrum and ∆ X,z i (cid:96) ( k ) = (cid:90) z i +∆ zz i − ∆ z dz dN X dz W ( z, z i )∆ X(cid:96) ( k, z ) , (2.8)where dN X dz is the source number density per redshift interval, W ( z, z i ) is a window functioncentered at z i with half-width ∆ z (with the integral of W ( z, z i ) dN X dz being normalized to– 3 –nity). Note that equation (2.7) follows the notation of reference [38], which reflects how thepublic code CLASS [39, 40] is built. The ∆ X(cid:96) ( k, z ) is the angular number count fluctuation ofthe X tracer, which is determined by density ( den ), velocity ( vel ), lensing ( len ) and gravity( gr ) effects [38, 41]: ∆ (cid:96) ( k, z ) = ∆ den (cid:96) ( k, z ) + ∆ vel (cid:96) ( k, z ) + ∆ len (cid:96) ( k, z ) + ∆ gr (cid:96) ( k, z ) . (2.9)The reader interested in the full expressions of the number counts fluctuations in equation(2.9) can find them in Appendix A. The relative importance between each of these termsdepends on the specific configuration (redshift bins, window functions, etc.) but some generalstatements can be made (see e.g., figure 1 of reference [40]). The main contribution is usuallygiven by the density term, whereas the gravity effects are subdominant (even by two ordersof magnitude at (cid:96) ∼ ), since it is relevant mostly at horizon scales. The lensing term isonly slightly scale-dependent, while the velocity one can be comparable to it at smaller scales,while stronger than it at lower multipoles (even almost one order of magnitude at very low (cid:96) ).This last statement holds especially for auto-bin correlations, while for power spectra amongdistant redshift bins the lensing term can overcome the velocity.The angular power spectra were computed using Multi_CLASS , the modified versionof
CLASS presented in [42, 43] which allows the user to compute cross-correlations betweendifferent tracers X (cid:54) = Y . We consider Gaussian window functions and fix the cosmologicalparameters from the Planck 2015 results with n s = 0 . , ln A s = 3 . , Ω cdm =0 . , Ω b = 0 . , h = 0 . [44].As can be seen by the above equations and Appendix A, there are four main ingredientsthat are needed to fully characterize one tracer X : • Redshift distribution dN X dz : the source number density per redshift interval is char-acterized by a shape which is a fundamental ingredient in the angular power spectracomputation (see equation (2.8)). Eventually, the number of sources in a specific red-shift bin is also necessary to compute the shot noise that enters in the estimate of theobserved ˜ C (cid:96) ’s in equation (2.3). More details about the redshift distributions of theconsidered tracers are given in sections 3.1.1 and 3.2.1. • Bias b X : it quantifies the mismatch between the distribution of matter and of the tracer X (see e.g., [45–52]). In this work we make use of the linear bias formulation. Indicatingthe local contrasts of matter and tracer X at position x respectively by δ ( x ) and δ X ( x ) ,we write δ X ( x ) ≡ n X ( x ) − ¯ n X ¯ n X = b X δ ( x ) , where n X is the comoving density of tracer Xand ¯ n X is its mean value. The bias appears as a linear factor in the density term ofequation (2.9) (see the full expression in appendix A). More details about the bias ofour tracers are given in sections 3.1.2 and 3.2.2. • Magnification bias s X ( z ) : it quantifies the change in the observed surface density ofsources of tracer X induced by gravitational lensing [53]. Two effects compete againsteach other: on one side the number of observed sources can increase due to a magni-fication of the received flux, which would make visible some sources right below thevisibility threshold (luminosity or magnitude for galaxies and Signal-to-Noise ratio forGWs); on the other side an increase of the area reduces the observed number density ofobjects. The magnification bias mainly affects the lensing term of equation (2.9), butenters also in the velocity and gravity terms. Further discussion about magnificationbias of our tracers can be found in sections 3.1.3 and 3.2.3.– 4 – Evolution bias f evo X : it reflects the fact that the number of elements of a tracer Xis not necessarily conserved in redshift due to the possible formation of new objects.The evolution bias can be written as [54–56]: f evog ( z ) = d ln (cid:18) a d N g dzd Ω (cid:19) d ln a , where a is thescale factor and d N X dzd Ω is the absolute distribution of objects of tracer X , which canusually be substituted by the observed distribution with good approximation [15]. Theevolution bias appears only in sub-leading contributions, since it is just present in thenon-dominant part of the velocity term and in the gravity term (which has a smallerinfluence with respect to the others) [40] of equation (2.9).We complete this section by explicitly summarizing the dependence of the various con-tributions to the angular number count fluctuations on the three biases types presented above(see appendix A for full expressions): ∆ den (cid:96) = ∆ den (cid:96) ( b X )∆ vel (cid:96) = ∆ vel (cid:96) ( s X , f evo X )∆ len (cid:96) = ∆ len (cid:96) ( s X )∆ gr (cid:96) = ∆ gr (cid:96) ( s X , f evo X ) (2.10)where dependencies on k and z are implied. In this section we describe how our two tracers (galaxies and GWs) are characterized, alongwith the theoretical frameworks used in their modeling. Since we deal with a statisticalapproach to count and describe the properties of galaxies and GWs, we use the notation dp/d Θ to express the probability distribution of a generic variable Θ . Our first tracers are actively star-forming galaxies. We do not deal with any specific galaxycatalog: we count and distribute galaxies on the basis of the observationally determinedSFRF at different redshifts (described in section 3.1.1). The SFR of a galaxy measures thestellar mass in solar units formed per year inside the galaxy. We briefly explain how it ismeasured in section 3.1.1. In this work we consider objects with three different SFR lowerlimits: ψ ≥ , , M (cid:12) / yr. The lower value of M (cid:12) / yr roughly corresponds to thelimit below which uncertainties in the star formation rate functions (SFRFs) are significant,especially at high redshift. The cut of M (cid:12) / yr is set to take into account the highly starforming dusty galaxies which constitute the bulk of the cosmic SFR. Finally, the highest limitof M (cid:12) / yr is set to take into account the most extreme star forming objects. The star formation rate functions SFRF = d N/d log ψ/dV correspond to the number den-sity of galaxies per cosmological comoving volume per logarithmic bin of SFR at a given cosmictime t or redshift z . In the last years several observations (e.g., UV+far-IR/submillimeter/radioluminosity functions and stellar/gas/dust mass functions) have allowed to robustly estimatethese functions. The SFR of a galaxy could in principle be estimated by its UV luminosity,since it is proportional to the quantity of young stars present in the galaxy. However, this– 5 –stimation can be easily biased by the presence of dust. Indeed, even a modest amount of dustcan significantly absorb the UV radiation and re-emit it in the far-IR/(sub)millimeter wave-lengths. Standard UV slope corrections (e.g., [57–59]) can still be applied to galaxies withrelatively low SFR ψ (cid:46) − M (cid:12) /yr , since the dust attenuation for them is mild. Therefore,deep UV surveys in the rest frame UV band are enough to robustly determine the SFRF atthe faint end. Instead, in highly star-forming galaxies with SFR ψ (cid:38) − M (cid:12) / yr dustobscuration is heavy and the corrections mentioned above are no more reliable (e.g., [60–64]).To soundly estimate their SFRF, it is necessary to use far-IR/(sub)millimeters observations.The latters have been exploited in many works over the recent years (e.g., [65–69]) to recon-struct, in combination with UV data, the SFRF for the whole SFR range at redshift z (cid:46) . Athigher redshifts, given the sensitivity limits of wide-area far-IR surveys, the reconstruction ofthe SFRF, especially at the bright end, is more uncertain. Useful information have been ob-tained from far-IR/(sub)millimeter stacking (see [70], [71]) and super-deblending techniques(see [72]), from targeted far-IR/(sub)millimeter observations (e.g., [73–75]) and from radiosurveys ([76]).All the above datasets have been fitted through simple Schechter functions by [77],obtaining: d Nd log ψ dV (log ψ, z ) = N ( z ) (cid:18) ψψ c ( z ) (cid:19) − α ( z ) e − ψ/ψ c ( z ) , (3.1)where the values of the redshift-dependent parameters N ( z ) , ψ c ( z ) and α ( z ) can be found intable 1 of [77] (see also figure 1 of [23]). From the SFRF we can obtain the number density ofgalaxies per unit comoving volume at different cosmic times t and the cosmic star formationrate density as: dNdV ( t ) = (cid:90) d log ψ d Nd log ψ dV (log ψ, t ) , (3.2) ρ ψ ( t ) = (cid:90) d log ψ ψ d Nd log ψ dV (log ψ, t ) . (3.3)Notice that the evolution of the cosmic star formation rate density with redshift, recon-structed in this way, is well in agreement with the available datasets (see [23]). The numberof galaxies per redshift bin dN/dz can be easily obtained multiplying equation 3.2 by thedifferential cosmological comoving volume dV /dz . The redshift distribution of our 3 galacticpopulations (galaxies with ψ > , , M (cid:12) / yr) can be obtained integrating the SFRFexcluding galaxies below a certain threshold ¯ ψ : dN ¯ ψ dz ( t, ψ ≥ ¯ ψ ) = dVdz (cid:90) ¯ ψ d log ψ d Nd log ψ dV . (3.4)In figure 1, left panel, the three dN ¯ ψ /dz of the galactic populations considered in thiswork are plotted as a function of redshift.From this quantity it is straightforward to compute the evolution bias for each galacticpopulation as f evog ( z ) = d ln (cid:16) a d N ¯ ψ ( z,ψ ≥ ¯ ψ ) dzd Ω (cid:17) d ln a . (3.5)– 6 – .1.2 Galaxy bias A fundamental ingredient entering in the computation of the observed number counts fluctu-ations in equation 2.9 is the bias b X ( z ) . Since our tracers are galaxies selected, counted anddivided by their SFR, we should connect the bias to this quantity. We adopt the procedure ofreference [78] associating the luminosity/SFR of the galaxy to the mass of the hosting darkmatter halo through an abundance matching technique and then assigning to a galaxy withgiven SFR the bias of the corresponding halo. Abundance matching is a standard method toderive a monotonic relationship between the galaxy and the halo properties by matching thecorresponding number densities in the following way: (cid:90) ∞ log ψ d log ψ (cid:48) d Nd log ψ (cid:48) dV = (cid:90) ∞−∞ d log M (cid:48) H d Nd log M (cid:48) H dV erf c (cid:26) log ( M H ( ψ ) /M (cid:48) H ) √ σ (cid:27) , (3.6)where M (cid:48) H is the halo mass, d N/d log M (cid:48) H /dV is the galaxy halo mass function i.e. themass function of halos hosting one individual galaxy (see Appendix A of [78]) and M H ( ψ ) isthe relation we are looking for. Finally, ˜ σ ≡ σ d log M H /d log ψ is the scatter around thatrelation (we set σ log ψ (cid:39) . following [79–81]). Once M H ( ψ ) is determined we assign toeach galaxy the bias corresponding to the halo associated to its SFR: b ( z, ψ ) = b ( z, M H ( z, ψ )) ,where b ( z, M H ) is computed as in [82] and approximated by [83]. Note that this formulationfor the bias is based on the excursion set approach. Other possible alternatives are present,such as the Effective Field Theory, Peak Theory, etc. (see e.g., [52, 84] and references therein).It is now easy to compute an effective bias for all the galaxies above a certain SFRthreshold ¯ ψ weighting b ( z, ψ ) by the corresponding galaxy distribution: b ¯ ψ ( z, ψ ≥ ¯ ψ ) = (cid:82) ∞ log ¯ ψ d log ψ d Nd log ψ dV b ( z, ψ ) (cid:82) ∞ log ¯ ψ d log ψ d Nd log ψ dV . (3.7)In figure 1, middle panel, we show the galaxy effective bias as a function of redshift forour galactic populations. We see, as expected, that it tends to increase with redshift and itis, in general, lower for galaxies with lower SFRs, since these are typically associated to lessmassive halos.Note that the bias computed in this section refers to star forming galaxies, which can bethe progenitors of quenched massive objects at the present time. We stress that, in order toestimate the bias for these galaxies, we should not look at their current SFR, but at the SFRof their progenitor when most of the stellar mass was accumulated, which is directly linkedto the mass of the host dark matter halo, via the abundance matching technique describedin equation (3.6). The magnification bias is another important factor entering in the angular power spectra com-putation. In fact, as we already described in section 2, the contribution to the angular numbercount fluctuations due to the lensing term of equation (2.9) can be comparable to that of thevelocity term (or even higher, especially when correlating objects between distant redshiftbins). As mentioned above, our galactic populations have SFR cuts ¯ ψ = 10 , , M (cid:12) / yr.The magnification bias for each of them is proportional to the logarithmic slope of their– 7 – N ψ /dz computed at ψ = ¯ ψ : s g, ¯ ψ ( z ) = − d log (cid:16) d N ψ ( z,>ψ ) dzd Ω (cid:17) d log ψ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ψ = ¯ ψ . (3.8)Using equation (3.4) we can show that the magnification bias can be directly related to theSFRF as: s g, ¯ ψ ( z ) = 25 ln 10 d Nd log ψ dV ( z, ¯ ψ ) dN ¯ ψ /dz dVdz . (3.9)Figure 1 (right panel) shows the magnification bias for our galactic populations as afunction of redshift. We can see that, at small z , the magnification bias decreases rapidly,especially for the tracers with higher SFRs: this is due to the fact that at small redshiftswe have less and less galaxies with high SFRs, therefore the function d N ( z, ψ ≥ ¯ ψ ) /dz/d Ω strongly depends on the choice of the faint end in SFR. Moreover, the overall magnificationbias for higher star forming galaxies is larger because they are less and a variation of the faintend SFR limit has a larger impact on their dN/dz . This is why, in general, the magnificationbias shape tends to be specular to the one of the dN/dz . z d N / d z ( z ) Redshift distribution z b ( z ) Bias ψ > M (cid:12) yr − ψ > M (cid:12) yr − ψ > M (cid:12) yr − z s ( z ) Magnification bias
Figure 1 . Full-sky redshift distributions ( left ), bias ( center ) and magnification bias ( right ) for allgalactic populations. Quantities referred to galaxies with Star Formation Rate ψ > M (cid:12) /yr, ψ > M (cid:12) /yr, ψ > M (cid:12) / yr are respectively in green, cyan and magenta lines. The other class of tracers that we are considering are GWs resolved signals coming from themerging of BH-BH, NS-NS and BH-NS binaries. In the next subsections, as already done forthe case of galaxies, we are going to illustrate how to compute their merging and detectedrates, their bias and magnification bias.
The merging rate of compact binaries is the number of merging per year per redshift bin d ˙ N /dz . Multiplying it by an observational time one gets the number of events per redshift The factor d ˙ N/dz is well defined since it stays put for typical observation times of the order of ∼
10 yr . – 8 –in in that period of time. Computing the merging and detected rates of compact binaries isa complex issue since it involves the necessity to understand and to correctly model many as-trophysical processes occurring on different time and spatial scales: from stellar astrophysics,to galaxy formation and evolution, to GW physics. However, a number of studies have ap-proached the problem combining population synthesis simulations (e.g., [85–91]) either withcosmological simulations (e.g., [92–96]) or with recipes on the cosmic star formation rate den-sity and metallicity distributions inferred from observations (e.g., [23, 97–102]). In this workwe follow the approach of reference [23], briefly sketched hereafter.The three main ingredients to compute the merging rates of compact binaries are: i) anobservational determination of the SFRF at different redshifts, ii) average chemical enrich-ment histories of individual galaxies and iii) outcomes from single stellar and binary evolutionsimulations.The first ingredient has already been described in section 3.1.1 and provides the galaxystatistics. In the following subsections we are going to describe respectively the other twoingredients and the way to combine them to compute the merging rates of compact binaries. Metallicity
The average chemical enrichment histories of individual galaxies is crucial since it allows toassociate a metallicity Z to galaxies with different properties (SFR, mass, age or morphologicaltype). Knowing the metallicity is fundamental because many binary evolutionary phenomenastrongly depend on it: stellar winds, supernova kicks, direct collapse, common envelope effects,etc. (for a more detail explanation of the main effects of metallicity on stellar and binaryevolution see section 5.1 and references therein). In [23] the chemical enrichment history of agalaxy with a given average SFR is reproduced with a simple model featuring a linear increaseof the metallicity in the early stages of the galaxy life up to a saturation value dependent onthe SFR. This model is an approximation of the more elaborated chemical evolution modelof [103] and [104] and well reproduces observations of both elliptical and disk galaxies (seee.g., [105–111]). The authors of [23] generate a metallicity distribution at given cosmic time t and SFR dp/d log Z (log Z | t, ψ ) taking into account the time spent by a galaxy with givenSFR in each bin of metallicity. Their final expression to compute the metallicity distributionis: dpd log Z (log Z | t, ψ ) = ∆ × ZZ sat ln(10)Θ H ( Z − Z sat ) + (1 − ∆) × δ D (log Z − log Z sat ) , (3.10)where Z sat ( t, ψ ) and ∆( t, ψ ) are parameters depending on the cosmic time and the SFR ofthe single galaxy. Z sat represents the saturation value of the metallicity and its typical valuesare in the range ∼ . − . Z (cid:12) , ∆ ∼ . − . specifies how quickly the metallicity saturates tosuch values as a consequence of the interplay between cooling, dilution and feedback processes. Stellar and binary evolution
The outcomes of single stellar and binary evolution simulations can provide three importantfactors: the remnant mass m • ( m (cid:63) , Z ) as a function of the zero age main sequence (ZAMS)star mass m (cid:63) and metallicity Z , a time delay distribution between the formation of the binaryand the merging dp/dt d and a mass ratio distribution dp/dq , where q is the ratio between theless and the more massive compact remnant.For the mass distribution of compact remnants we take as a reference the m • ( m (cid:63) , Z ) relation given in [89], and we generate a probability distribution function just applying a– 9 –ogarithmic gaussian scatter ( σ = 0 . dex) around the m • ( m (cid:63) , Z ) value to take into accountpossible uncertainties coming from stellar evolutionary processes as in [23]: dpd log m • ( m • | m (cid:63) , Z ) = 1 √ πσ exp (cid:2) − (log m • − log m • ( m (cid:63) , Z )) / σ (cid:3) . (3.11)To get rid of the dependence on the initial stellar mass it is sufficient to select an initial massfunction (IMF) φ ( m (cid:63) ) (in this work we used the Chabrier one [112] ) and integrate over theinitial stellar masses weighting the integral with the IMF: dpd log m • ( m • | Z ) = (cid:90) ¯ m (cid:63) dm (cid:63) φ ( m (cid:63) ) dpd log m • ( m • | m (cid:63) , Z ) , (3.12)where ¯ m (cid:63) ∼ M (cid:12) is the ZAMS star mass limit originating a NS remnant . However, sincethe amplitude of the gravitational wave events is determined by the chirp mass M •• ≡ m • q / / (1 + q ) / , rather than by the primary mass m • , we should make use of the massratio distribution to change variable and determine the probability distribution function fora given chirp mass in the following way: dpd M •• ( M •• | Z ) = (cid:90) dq dpdq dpdm • ( m • ( M •• , q ) | Z ) dm • d M •• ( q ) , (3.13)where m • ( M •• , q ) = M •• (1 + q ) / /q / , dm • /d M •• = (1 + q ) / /q / and the distribution dp/dq is taken from binary evolution simulations. In particular, the mass ratio distribution forBH-BH mergers scales linearly with q : dp/dq ∝ q (see e.g., [87, 97, 113]), instead, for NS-NSor BH-NS merging, the mass ratio distribution tends to be flatter (see [86, 95, 114–116]). Weset the q distributions ranges for the three types of merging events on the basis of the allowedmasses for each CO type. The authors of reference [23] have checked that the merging ratedepends very little on the chosen q distribution.The last ingredient provided by observations and simulations is the probability distri-bution function for the time delay between the formation of the binary and its merging: dp/dt d ∝ t − d [91, 114], normalized to unity between a minimum value of t d,min ∼ Myr andthe age of the universe.
Computing merging and detected rates
After this brief overview of the main ingredients, we can compute the cosmic merging ratesof compact remnant binaries per redshift and chirp mass interval in the following way: d ˙ N merge dz d M •• ( t, M •• ) = f eff dVdz (1 + z ) (cid:90) dt d dpdt d (cid:90) d log ψ d N (log ψ, t − t d ) d log ψ dV ψ ×× (cid:90) d log Z dpd log Z (log Z | t − t d , ψ ) dpd M •• ( M •• | Z ) . (3.14) In principle the IMF could be different and also dependent on the galaxy properties as SFR or metallicity.However the authors of reference [23] showed that its choice has only a mild impact on the detected GWs rate. Note that in equation (3.12), the integral should contain the quantity φ ( m (cid:63) ) / (cid:82) dm (cid:63) φ ( m (cid:63) ) m (cid:63) . How-ever, in literature the denominator is usually left implicit because of the IMF normalization condition (cid:82) dm (cid:63) φ ( m (cid:63) ) m (cid:63) = 1 M (cid:12) , though the reader should keep track of the measurement units (i.e. the factorcomputed in equation (3.12) is a distribution of remnant masses per unit of star formed mass). The primary mass m • is the mass of the most massive compact remnant of the binary. – 10 –o grasp the meaning of this cumbersome expression let us look first at the innermostintegral: it represents the probability of formation of a compact remnant with chirp mass M •• in a galaxy with SFR ψ . This quantity is then integrated over all the galaxies weightedwith the galaxy statistics (the SFRF). Since we are computing the merging rates at the cosmictime t , the quantities related to the formation of the binary should be computed at t − t d , with t d being the time delay between the formation and the merging of the binary. The outermostintegral over the time delay is finally performed to account for all the time delays. The factor dV /dz is just the differential comoving volume, while the (1 + z ) factor in the denominatorkeeps into account the cosmological time dilation. Finally, the factor f eff is defined as thefraction of primary compact remnants hosted in binary systems with characteristics apt toallow merging within a Hubble time.The factor f eff is the result of many different and complex physical processes relatedto stellar and dynamical evolution (binary fraction, common envelope development/survival,SN kicks, mass transfers, etc.), so it could in principle depend on metallicity and binarytype. In stellar and binary evolution simulations (e.g., [86, 89, 91, 97, 102, 117, 118]) thisquantity is naturally obtained, but largely dependent on model assumptions. This is why,at this stage, as already done in other works (e.g. [99], [101], [19]), we set it empiricallyby normalizing the local BH-BH, BH-NS and NS-NS merger rates to the logarithmic averagevalues of the confidence interval measured by the LIGO/Virgo collaboration after the O1and O2 runs (see [119]):
30 Gpc − / yr for BH-BH,
650 Gpc − / yr for NS-NS and
25 Gpc − / yrfor BH-NS (this last choice is less certain since the BH-NS local merging rate is limited onlyby an upper value). We stress that in doing so the factor f eff loses all the possible metallicitydependence, which, however, is highly uncertain and model dependent (we will come back toits metallicity behaviour in section 5). Therefore, this factor acts only on the normalizationof the merging rates and can be changed when a more accurate determination of the localrates will be done after further GWs observations. Thus, the difference in the merging ratesnormalization between the three types of merging events is given by the different f eff factors,while the difference in shape is given by the distribution dp/d M •• ( M •• | Z ) in equation (3.13),depending on the stellar and binary evolution prescriptions.Once the merging rates per chirp mass bin are computed, it is easy to derive the detectedrates by a specific GW detector. As mentioned in section 1, we consider the ET instrument.The rates per unit redshift, chirp mass and signal to noise ratio (SNR) can be computed as: d ˙ N merge dz d M •• dρ ( ρ | z, M •• ) = d ˙ N merge dz d M •• dpdρ ( ρ | z, M •• ) (3.15)where dp/dρ is the probability distribution of SNR dependent on redshift, chirp mass and onthe sensitivity curve of the detector (for a full treatment of the dp/dρ see [23, 101, 120, 121]).Therefore, the rates with SNR ρ > ¯ ρ can be computed as: d ˙ N ¯ ρ dz ( z, ρ ≥ ¯ ρ ) = (cid:90) ¯ ρ dρ d ˙ Ndz dρ ( ρ, z ) = (cid:90) ¯ ρ dρ (cid:90) d M •• d ˙ Ndz d M •• dρ ( M •• | z, ρ ) . (3.16)We consider the GW event detected when the SNR is higher than ¯ ρ = 8 . The detectedrates by ET for BH-BH, NS-NS and BH-NS are shown in figure 2 (left panel).Eventually, the rate for GWs events above a given SNR can be used to compute theevolution bias: f evoGW ( z ) = d ln (cid:16) a d ˙ N ¯ ρ ( z,ρ ≥ ¯ ρ ) dzd Ω (cid:17) d ln a . (3.17)– 11 – .2.2 Bias for GW events Since we consider GWs produced by the merging of COs of stellar origin, their signals origi-nate from galaxies and trace their distribution, so that they trace the underlying total matterdistribution the same way their host galaxies do. For this reason, GWs events can be charac-terized by the same bias of their hosts. Since the COs mergers take place in all galaxy typeswith different rates, the correct estimate of their bias needs to take into account which galaxytypes are contributing most/less to the detected mergers, giving proportioned weight to theirbias values when estimating that of all GWs events.In order to assign a redshift dependent bias to the GW events, we make use of the bias b ( z, ψ ) , computed in section 3.1.2, associated to a galaxy at a given redshift with given SFRand we weight it through the quantity d ˙ N merge /dz/dρ/d log ψ which keeps into account thecontribution of the different SFRs (i.e. of different galaxies) to the total merging rates at agiven redshift and SNR. This differential merging rate can be computed from equation (3.14)and (3.15) not integrating over the SFR. Therefore, to compute the bias for gravitationalwaves we use the following expression: b GW ( z, ρ ) = (cid:82) d log ψ d ˙ Ndz dρ d log ψ b ( z, ψ ) (cid:82) d log ψ d ˙ Ndz dρ d log ψ . (3.18)The effective bias, i.e. the bias for GWs with a SNR above a certain threshold ¯ ρ , is noweasy to compute: b GW, ¯ ρ ( z, > ¯ ρ ) = (cid:82) ¯ ρ dρ d ˙ Ndz dρ b ( z, ρ ) (cid:82) ¯ ρ dρ d ˙ Ndz dρ . (3.19)The bias for the detected events ( ¯ ρ = 8 ) is shown in figure 2 (middle panel). Theinterpretation of the shape of the GW bias is not trivial and explained in the following. Atlow redshift its value is ∼ since the only galaxies that contribute to the GW signals havelow SFR and consequently a smaller bias. The following rapid increase with redshift is dueto two factors: the first is just the standard growth with redshift of the galaxy bias, thesecond is that, increasing the redshift, there are more and more highly star forming galaxiesthat contribute to the GW events. These galaxies, as shown in figure 1, are more biased.At redshift z ∼ the GW bias flattens. Again, this is due to different astrophysical effects.In particular, the redshift increase of the galaxies bias is compensated by the fact that athigh redshift the detected GW events receive a larger contribution by less star forming and,thus, less biased galaxies. This is due to two facts: firstly, the number of highly star forminggalaxies tends to decrease at redshift z (cid:38) − ; secondly, in galaxies with high SFR themetallicity is also high and, consequently, the compact remnants produced are less massive.This means that galaxies with larger SFRs tend to produce GW events with lower chirp mass(see section 5.1). However, at high redshift the detector starts not to see anymore these lowchirp mass events and, due to this selection effect, the GW events detected at higher andhigher redshifts come from galaxies with lower SFR and are, consequently, less biased.– 12 – .2.3 Magnification bias for GW events Similarly to the galaxy case and as done in ref. [15], the magnification bias for GW eventswith ρ > ¯ ρ is the logarithmic slope of their dN ρ /dz ( z, > ρ ) computed at ρ = ¯ ρ : s GW, ¯ ρ ( z ) = − d log (cid:16) d ˙ N ρ ( z,>ρ ) dz d Ω (cid:17) dρ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ =¯ ρ , (3.20)which, after some algebraic manipulation, can be rewritten as: s GW, ¯ ρ ( z ) = ¯ ρ d ˙ N ( z, ¯ ρ ) dz dρ d ˙ N ¯ ρ /dz . (3.21)We shown in figure 2 (right panel) the magnification bias for detected mergers ( ¯ ρ = 8 ). Itcan be seen that the magnification bias for NS-NS events features a fast growth with redshiftbecause the NS-NS distribution in SNR is peaked at lower values of ρ with respect to BH-BHor BH-NS events (see figure 3). So, as the redshift increases, the peak of such distributionshifts toward values ρ (cid:46) : the choice of the faint end of SNR has then a huge effect onNS-NS events. Instead, for BH-BH and BH-NS events, the distribution in SNR ratio is muchbroader, even at high redshifts: the choice of the faint end of SNR has not a large impacton the number of detections. For this reason the magnification bias for those events alwaysremains at moderate values.In figure 3 we show the SNR probability distribution functions for BH-BH, BH-NS andNS-NS events at redshift z = 0 . left panel, z = 1 middle panel and z = 2 right panel. z − d N / d z ( z ) Redshift distribution z b ( z ) Bias
BH-BH BH-NS NS-NS z s ( z ) Magnification bias
Figure 2 . Full-sky redshift distributions for an observation time T obs = 1 yr ( left ), bias ( center )and magnification bias ( right ) for all GWs tracers, as detected by ET. Quantities referred to BH-BH,BH-NS, NS-NS mergers are respectively in red, blue and yellow lines. A useful thing to notice is that changing the cosmological parameters values affectsthe description of the tracers (both galaxies and GWs) only as a volume term dV /dz in thecomputation of the redshift distributions, as can be seen in equations (3.4) (for galaxies)and (3.14) (for GWs). This implies that, changing the cosmological parameter values, all theredshift distributions will vary in the same way, making it difficult to constrain them.– 13 –
10 20 300 . . . d p / d ρ z = 0 . ρz = 1 . BH-BH BH-NS NS-NS z = 2 . Figure 3 . SNR ρ normalized distributions at redshifts z = 0 . ( left ), z = 1 . ( center ) and z = 2 . ( right ) for all GWs tracers. Distributions for BH-BH, BH-NS, NS-NS are respectively in red, blueand yellow lines. The black dashed vertical lines correspond to the limit of ρ = 8 . Making use of the formalism presented in section 2 we compute the cross-correlation angularpower spectra C (cid:96) ’s between the tracers presented in section 3. Note that the GWs eventsconsidered in the cross-correlation are computed from the whole galaxy distribution, not justfrom galaxies with the specified SFR cuts, in order to realistically take into account all theGWs signals that can be detected by a specific instrument (ET in our case). We consideredboth tomographic and non tomographic approaches, whose details and results are providedin sections 4.1 and 4.2 respectively.Following ref. [15], we can organize the angular power spectra from different tracers andredshift bin couples in a data vector C (cid:96) ordered as C (cid:96) = C gg (cid:96) ( z , z ) ... C gGW (cid:96) ( z , z ) ... C GWGW (cid:96) ( z , z ) ... (4.1)where g and GW respectively refer to our galaxy and gravitational wave tracers. Given N bin redshift bins, the C (cid:96) is a N dimensional data vector, where its I th element can be associatedto two indices ( I , I ) , corresponding to the two tracers and redshift bins of the angular powerspectra in that specific entry. As an example, the first ( I = 1 ) entry is associated to the coupleof indices [ I = g z , I = g z ] . From this, one can write the covariance matrix (Cov (cid:96) ) IJ , whoseelements are given by (Cov (cid:96) ) IJ = ˜ C I J (cid:96) ˜ C I J (cid:96) + ˜ C I J (cid:96) ˜ C I J (cid:96) , (4.2)where the ˜ C (cid:96) are the angular power spectra of equation (2.3).In order to characterize the magnitude of the signal that could be extracted by thecross-correlations and determine whether it could be discerned from the noise, it is useful tocompute a Signal-to-Noise ratio (S/N). With this purpose, we compute two types of S/N.– 14 –he first one is an estimate of the S/N of the C (cid:96) ’s for each combination of redshift bins(at fixed tracers couple). This provides a total of N bin × N bin S/N values. While on the onehand they do not take into account correlations with the other redshift bins (which can benon negligible, especially due to lensing effects), on the other hand they provide an unpackedinformation about which bin combinations are the most powerful in terms of signal. At fixedmultipole (cid:96) , the S/N computed in this way can be written as: (cid:18) SN (cid:19) I ,I ] ( (cid:96) ) = f sky (2 (cid:96) + 1) (cid:0) ˜ C [ I ,I ] (cid:96) (cid:1) σ I ,I ] ( (cid:96) ) = f sky (2 (cid:96) + 1) (cid:0) ˜ C [ I ,I ] (cid:96) (cid:1) (cid:34) ˜ C [ I ,I ] (cid:96) ˜ C [ I ,I ] (cid:96) + (cid:0) ˜ C [ I ,I ] (cid:96) (cid:1) (cid:35) (4.3)where σ I ,I ] ( (cid:96) ) is obtained imposing I = J in equation (4.2).The second method for the computation of the S/N provides (still at fixed tracers couple)one single S/N estimate for the whole probe, taking into account also the covariance betweenthe C (cid:96) ’s of different redshift bins. It is computed as: (cid:18) SN (cid:19) ( (cid:96) ) = f sky (2 (cid:96) + 1) · C T (cid:96) · Cov − (cid:96) · C (cid:96) (4.4)Note that, in order to compute the S/N of the cross-correlation, we do not make use ofthe auto-correlation C (cid:96) ’s appearing in the vector of equation (4.1). In this way we compute theS/N related only to the cross-correlation part, avoiding the contribution from auto-correlationswhich is likely to increase the S/N, if the auto-correlations are free of systematic effects. Therationale behind this treatment is to be conservative and assume that the cross-correlationsignal is less prone to systematic effects compared to the auto-correlations, as it is indeed thecase for other Large Scale Structure tracers.More often, a cumulative Signal-to-Noise ratio (cid:32) SN (cid:33) ( (cid:96) < (cid:96) max ) is considered. In bothcases it is defined as (cid:18) SN (cid:19) ( (cid:96) < (cid:96) max ) = (cid:118)(cid:117)(cid:117)(cid:116) (cid:96) max (cid:88) (cid:96) (cid:48) = (cid:96) min (cid:18) SN (cid:19) ( (cid:96) (cid:48) ) . (4.5)In sections 4.1 and 4.2 we show the results for the two S/N calculations. Note that in thenon tomographic case, characterized by one single redshift bin, the two estimates coincide. We cross-correlate all the considered galaxy tracers (i.e. galaxies with ψ > , , M (cid:12) / yr)with all the GWs tracers (i.e. BH-BH, BH-NS, NS-NS) along three or four redshift bins z i,j .The number of bins and their ranges differ from case to case due to the different redshiftranges in which these tracers can be defined (see section 3). In table 1 we provide the redshiftbinning considered for each probe. In this section we provide results for the exemplificativecase of f sky = 0 . and up to a maximum multipole of (cid:96) max = 100 , corresponding to the bestangular resolution reachable by ET (see e.g., [122]). The angular power spectra for all thetracers combinations are shown in Appendix B.In figures 4, 5, 6 we provide the cumulative S/N computed from the first method (i.e.applying equation (4.5) to (4.3)) for all combinations of tracers and redshift bins. Eachsubplot shows the estimates for one specific galactic tracer, with every GW tracer. The– 15 –WBH-BH BH-NS NS-NS ψ > M (cid:12) / yr ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ g ψ > M (cid:12) / yr . ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ . ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ . ≤ z ≤ ≤ z ≤ ≤ z ≤ ψ > M (cid:12) / yr ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ ≤ z ≤ Table 1 . Redshift binning for the tomographic case. following notation regarding the redshift binning is adopted: given a g × GW couple, thenotation z i − z j means that we are cross-correlating galaxies in bin z i with GWs in bin z j .Note that each subplot refers to cross-correlations between redshift bins z i − z j , but each z i,j can actually be different between different tracers according to table 1. For this reason,curves with different colors should not be directly compared to one another, since they referto different ranges.As for the S/N values, it can be seen that in several cases S / N( < (cid:96) max ) > . In particular: i) for all galaxy tracers, the highest S / N( < (cid:96) max ) are found for correlations among the sameredshift bins; ii) cross-correlations between distant redshift bins are those providing a lower S / N( < (cid:96) max ) : even though effects such as lensing can induce even a strong correlation betweendistant objects, in many of the cases considered here it is not enough to strengthen the S / N ; iii) considering correlations among the same bins, the S / N( < (cid:96) max ) are strongly sensitive tothe amount of detected sources: for large redshift values (bins z , , ) they are always higherin the BH-BH case, followed by BH-NS and eventually by NS-NS. Indeed, the BH-BH casecorresponds to a higher number of merging events (as can be seen by looking at its redshiftdistribution of the left panel of figure 2). This provides a smaller amount of shot noise, whichcontributes to making the S / N of equation (4.3) higher. On the other hand, at low redshift(bin z ) the redshift distribution of NS-NS mergers is significantly higher than the others,which is reflected in a S / N( < (cid:96) max ) which is often bigger; iv) in analogy with the previouspoint, at fixed GW tracer the higher is the cut in SFR, the lower is the S/N, reflecting thesmaller number of galaxies considered. – 16 – z − z z − z BH-BH BH-NS NS-NS z − z z − z z − z z − z z − z z − z S / N ( < ‘ m a x ) z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z ‘ max z − z z − z z − z z − z z − z z − z galaxies: ψ > M (cid:12) / yr Figure 4 . Cumulative Signal-to-Noises S / N( < (cid:96) max ) (equations (4.3) and (4.5)) for the cross-correlations cases between galaxies with ψ > M (cid:12) / yr with all three types of GWs signals (BH-BHin red, BH-NS in blue, NS-NS in yellow). Horizontal dashed lines correspond to S / N( < (cid:96) max ) = 1 .The plot refers to T obs = 1 yr and f sky = 0 . . – 17 – z − z z − z BH-BH BH-NS NS-NS z − z z − z z − z z − z z − z z − z S / N ( < ‘ m a x ) z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z z − z ‘ max z − z z − z z − z z − z z − z z − z galaxies: ψ > M (cid:12) / yr Figure 5 . Cumulative Signal-to-Noises S / N( < (cid:96) max ) (equations (4.3) and (4.5)) for the cross-correlations cases between galaxies with ψ > M (cid:12) / yr with all three types of GWs signals (BH-BHin red, BH-NS in blue, NS-NS in yellow). Horizontal dashed lines correspond to S / N( < (cid:96) max ) = 1 .The plot refers to T obs = 1 yr and f sky = 0 . . – 18 – z − z z − z BH-BH BH-NS NS-NS z − z z − z z − z z − z S / N ( < ‘ m a x ) z − z z − z z − z z − z z − z z − z z − z z − z ‘ max z − z z − z z − z z − z galaxies: ψ > M (cid:12) / yr Figure 6 . Cumulative Signal-to-Noises S / N( < (cid:96) max ) (equations (4.3) and (4.5)) for the cross-correlations cases between galaxies with ψ > M (cid:12) / yr with all three types of GWs signals (BH-BHin red, BH-NS in blue, NS-NS in yellow). Horizontal dashed lines correspond to S / N( < (cid:96) max ) = 1 .The plot refers to T obs = 1 yr and f sky = 0 . . The cumulative S/N for the whole probe, computed applying equation (4.5) to (4.4), isshown in figure 7. The line-styles refer to a specific galaxy tracer, while the colors distinguishbetween GW types (as indicated in the legend). It can be seen that generally it reachesvalues above unity for all the tracers combinations. Note that the cross-correlation signalovercomes the noise already for relatively low multipoles, at around (cid:96) max ∼ − . The S/Nis particularly high especially for the ψ > M (cid:12) / yr × NS-NS case, where it reaches a value of ∼ . The NS-NS case is also the most dependent on the chosen SFR cut, because the peakof the detected NS-NS distribution is at rather low redshift ( z ≤ ), where the distributionof highly star forming galaxies tends to fall down (see figures 1 and 2, left panels). We stress– 19 –WBH-BH BH-NS NS-NS ψ > M (cid:12) / yr ≤ z ≤ ≤ z ≤ ≤ z ≤ g ψ > M (cid:12) / yr . ≤ z ≤ . ≤ z ≤ . ≤ z ≤ ψ > M (cid:12) / yr ≤ z ≤ ≤ z ≤ ≤ z ≤ Table 2 . Redshift binning for the non-tomographic case. again that the different cases should not be directly compared since they refer to differentredshift ranges, depending on the tracers characteristic intervals. All in all, cross-correlationsbetween the treated tracers, adopting a tomographic approach, can be informative given therather high S/N ratio values. ‘ max S / N ( < ‘ m a x ) ψ > M (cid:12) / yr × BH-BH ψ > M (cid:12) / yr × BH-NS ψ > M (cid:12) / yr × NS-NS ψ > M (cid:12) / yr × BH-BH ψ > M (cid:12) / yr × BH-NS ψ > M (cid:12) / yr × NS-NS ψ > M (cid:12) / yr × BH-BH ψ > M (cid:12) / yr × BH-NS ψ > M (cid:12) / yr × NS-NS
Figure 7 . Cumulative Signal-to-Noises S / N( < (cid:96) max ) for all cross-correlation cases (equations (4.4)and (4.5)). The horizontal dashed line corresponds to S / N( < (cid:96) max ) = 1 . Line-styles refer to galaxies( ψ > M (cid:12) / yr in full line, ψ > M (cid:12) / yr in dashed line, ψ > M (cid:12) / yr in dotted-dashed line) whilecolors refer to gravitational waves (BH-BH in red, BH-NS in blue, NS-NS in yellow). The plot refersto T obs = 1 yr and f sky = 0 . . In this subsection we compute the cross-correlations between our tracers without a tomo-graphic approach. This is done to see how the measured cross-correlation signal can differwhen squeezing all the detected sources into one single bin, without a sliced tomographicapproach. As for the previous case, the redshift ranges of all the tracers combinations differfor each combination. In table 2 we provide the redshift ranges considered for each probe.In figure 8 (left panel) we show the angular power spectra for all combinations of tracers.As in figure 7, line-styles/colors refer to galaxies/GWs. It can be seen that the power spectraare higher when galaxies with a higher SFR are considered, as logically expected, since a– 20 –igher SFR leads to a larger absolute number of remnants, and so to a larger amount ofmerging pairs.In figure 8 (right panel) we provide the cumulative S / N( < (cid:96) max ) for each probe. Westress again that the S/N, computed applying equation (4.5) to equations (4.3) and (4.4), arein this case coincident. It can be seen that for an (cid:96) max large enough a S / N( < (cid:96) max ) > isalways reached for any tracer combination. The cross-correlation signal overcomes the noiseat around (cid:96) max ∼ − . The NS-NS contribution in this case is lower with respect to thetomographic approach because, even if their shot noise is small, their C (cid:96) ’s values are alsosmall, due to the fact that NS-NS mergers are mostly seen at low redshifts making theirdistribution rather different with respect to the galaxies one.By comparing figure 7 and figure 8 (right panel) it is possible to gauge the fact thatassuming a tomographic approach is indeed an advantage, since the extra (radial) informationprovided contributes to build a stronger S/N. ‘ . . . . . C ‘ × − ‘ max . . . . . . . . S / N ( < ‘ m a x ) ψ > M (cid:12) / yr × BH-BH ψ > M (cid:12) / yr × BH-NS ψ > M (cid:12) / yr × NS-NS ψ > M (cid:12) / yr × BH-BH ψ > M (cid:12) / yr × BH-NS ψ > M (cid:12) / yr × NS-NS ψ > M (cid:12) / yr × BH-BH ψ > M (cid:12) / yr × BH-NS ψ > M (cid:12) / yr × NS-NS
Figure 8 . Left: angular power spectra C (cid:96) ’s for all cross-correlations cases. Right: cumulative Signal-to-Noises S / N( < (cid:96) max ) for all cross-correlations cases. The horizontal dashed line corresponds to S / N( < (cid:96) max ) = 1 . Line-styles refer to galaxies ( ψ > M (cid:12) / yr in full line, ψ > M (cid:12) / yr in dashedline, ψ > M (cid:12) / yr in dotted-dashed line) while colors refer to gravitational waves (BH-BH in red,BH-NS in blue, NS-NS in yellow). The plot refers to T obs = 1 yr and f sky = 0 . . In this section we are going to discuss about the possibility to exploit GW × LSS cross-correlations to compare and test the validity of different astrophysical scenarios concerningthe formation, evolution and merging of COs binaries. Given the uncertainties in the astro-physics and the enormous modeling possibilities, it is at the moment unlikely to be able tounequivocally determine the validity of one specific combination of prescriptions with respectto any possible other. For this reason, our approach does not aim at excluding or validatingspecific models, whereas it looks at the possibility to apply this methodology to some gen-eral proof-of-concept cases, leaving the chance to use this formalism to whom is interestedin applying it for specific tests. In fact, it is worth stressing again that the application we– 21 –resented in this section is not the only possible one. Any other astrophysical formulation andmodeling that influences the estimate of redshift distributions and/or biases of the consideredtracers can in principle be addressed. This can cover a wide range of possibilities, from thestellar modeling (especially regarding the estimation of the CO mass and the constrain of allthose processes influencing it) to galaxy evolution and SFRs calculation methods.To compare two models (which we describe in section 5.1) and investigate the possibilityto discern them, we make use of a ∆ χ statistics, whose formalism and results are presentedin section 5.2. The two cases that we compare in this study differ in the treatment of the effect of thegalaxy metallicity on binary evolution. The metallicity value of the ZAMS stars constitutingthe binary, as already mentioned in section 3.2.1, is fundamental in order to determine thesubsequent binary evolution. The value of metallicty has a strong impact both on the COsmass and on the number of COs merging per unit of star forming mass. In particular thehigher is the metallicity the lower is the chirp mass and the number of merging binaries perunit of star forming mass. The physical reasons of these dependencies are well explained ine.g., [97], [114], [118].In section 3.2.1 and throughout this paper we accounted for the differences in the com-pact remnant masses due to the different metallicities present in galaxies through the factor dp/d log M •• ( M •• | Z ) in equation (3.14). However, as explained in section 3.2.1, we didnot consider the possible dependence on metallicity of the merging efficiency due to binaryevolution effects. This is translated in the fact that we have chosen a factor f eff independenton metallicity, whose value was determined by normalizing the merging rates to the localvalue constrained by the LIGO/Virgo team.In this specific section we change approach, aiming to address the issue of the possibilityto distinguish, through GW × LSS cross-correlations, between different astrophysical modelsof stellar and binary evolution. Therefore, we compare our precedent results obtained with aconstant f eff (we refer to it as "benchmark case") with one of the models with a metallicitydependent f eff ( Z ) . As an example we choose the reference model presented in [118]. Themetallicity dependence of the number of merging events per unit star forming mass of such amodel is shown in figure 1 of [116], thin lines (we refer to it as "Z-dependent case"). However,the reader should keep in mind that we are considering only the shape of this factor as afunction of metallicity, the normalization is still fixed by the local merging rate given by theLIGO/Virgo teams. The equation to compute the merging rates now becomes: d ˙ N merge dz d M •• ( t, M •• ) =˜ a dVdz (1 + z ) (cid:90) dt d dpdt d (cid:90) log ψ d N (log ψ | t − t d ) d log ψ dV ψ ×× (cid:90) d log Z dpd log Z (log Z | t − t d , ψ ) f eff ( Z ) dpd M •• ( M •• | Z ) , (5.1)which is identical to equation (3.14), except for the metallicity dependent f eff factor, whichnow takes part in the integration. The factor ˜ a instead guarantees the normalization to theLIGO/Virgo local values mentioned above.Thus, we are comparing a simple case (the benchmark one), where the dependence onmetallicity enters only in the remnant mass distribution, but not in the number of mergingbinaries per unit of star forming mass, with a more realistic model (the Z-dependent case)– 22 –here this latter dependence is included. The differences in the shape of the merging ratesbetween the two cases enters in the computation of the C (cid:96) ’s allowing the possibility to distin-guish between them. On the other hand differences in the absolute number of sources affectsthe shot noise of the two cases, in particular the number of merging events of the benchmarkcase tends to be, on average, higher with respect to the Z-dependent case, resulting in a lowershot noise.Clearly this is a case study to check whether it is possible to detect this metallicityimprinting through GW × LSS cross correlations, but this technique could be in principlepursued in more refined studies to test different astrophysical models. One of the reasons ofexploiting the cross-correlation formalism to test these two cases is given by the fact that,since SFR and metallicity are interconnected parameters, a dependence on the metallicityof the efficiency with which COs binaries merge will (non-trivially) be expressed also as adependence on the SFR, causing GWs mergers to correlate differently with galaxies of differentSFR values. We finally remark that we are going to compare these two models only for theBH-BH merging case for two main reasons: the first one is that the metallicity dependenceof the number of merging events is stronger for the BH-BH case (see [116]), the second one isthat BH-BH events are much more frequent with respect to the other types of merging: this,as already seen, reduces the shot noise enhancing the S/N of the cross-correlation.
In this section we provide the forecasts for discerning the benchmark scenario from the metal-licity dependent one. We make use of a ∆ χ statistics to evaluate a S/N, whose value(above/below unity) can provide information on how different the two models (one called asFiducial and the other as Alternative) are. Following the same approach of [15] we define aS/N as: (cid:18) SN (cid:19) ∼ ∆ χ := f sky (cid:96) max (cid:88) (2 (cid:96) + 1)( C Alternative (cid:96) − C Fiducial (cid:96) ) T Cov − (cid:96) ( C Alternative (cid:96) − C Fiducial (cid:96) ) , (5.2)where C Fiducial / Alternative (cid:96) is a vector containing the C (cid:96) ’S from the Fiducial/Alternative model,organized with the same logic of equation (4.1) and where the Cov (cid:96) is the covariance matrix asin equation (4.2), built with the C (cid:96) ’s of the fiducial model. Since the entries of the covariancematrix depend on which model is assumed as fiducial, the final forecasts also depend on thischoice. For this reason, we computed S/N in both cases and compared them.In figure 9 we provide the S/N obtained considering galaxies with ψ > , , M (cid:12) / yr.We show results for different observed sky-fractions f sky , observation times T obs and for bothmodels assumed as fiducial. More precisely, our results are expressed not only as a functionof T obs (on the horizontal axis) but, instead, of the product r · T obs . The quantity r is amultiplicative fudge factor to the merging rate of GWs introduced to take into account anypossible uncertainty in the modeling of this quantity. Note that r and T obs are degenerate:for example, observing for T obs = 1 yr with a factor r = 2 yields the same result as observingfor T obs = 2 yr with a factor r = 1 . The case of r = 1 corresponds to the scenario in whichthe models used here are the "true" ones. It is worth noticing that the r factor has the sameeffect of the ˜ a factor that quantifies the normalization to the local observed rate introduced inequation 5.1, since they both are multiplicative factors to the merger rate. For this reason, the r factor can also be seen as absorbing the uncertainties on the local merging rates estimates.– 23 –irst of all, it can be seen that when the benchmark model is assumed as fiducial, theforecasts are significantly better compared to the opposite (for fixed f sky and T obs ): this isdue to the fact that this model predicts a higher number of GWs mergers, providing a smallershot noise contribution. For analogous reasons, when comparing the panels in figure 9, itcan be seen that results are more optimistic when considering galaxies with ψ > M (cid:12) / yr:the higher the number of sources (galaxies in this case) the better the results. This is alsoreflected on the fact that the case with galaxies of ψ > M (cid:12) / yr is the most pessimistic.Looking in detail at each of the figures, we can see that in the ψ > M (cid:12) / yr case a S/Nabove unity can be reached in a relatively short time: even for small observed fractions of thesky (e.g., f sky = 0 . ) not more than 3 years of observation would be required to marginallydistinguish the two scenarios. Looking instead at the most pessimistic case, in which onlygalaxies with ψ > M (cid:12) / yr are considered, approximately 5 to 10 years of observation wouldbe required to reach S / N ∼ in the case of benchmark model assumed as fiducial. A muchhigher observation time (at least above 10 years) is required when assuming the Z -dependentscenario as fiducial. The case of ψ > M (cid:12) / yr lies in between, with a still fairly optimisticprediction.All in all our results are rather promising, especially when considering galaxies with ψ >
10 M (cid:12) / yr and ψ >
100 M (cid:12) / yr : if the benchmark case is the fiducial one, deviations fromit can be detected after just of observational time. If the Z -dependent case is the fiducial,we should be able to detect variations from it in (cid:46) of observations. For highly starforming galaxies with ψ >
300 M (cid:12) / yr instead, some more time is required to distinguish thetwo models. Still an observational time (cid:46)
10 yr is enough if the benchmark case is consideredas fiducial.Finally, we stress again that this forecast does not aim at testing or excluding anyof the two models considered here. It aims instead at showing how different astrophysicalprescriptions (such as a Z dependence on the f eff factor) could in principle be distinguishedthrough the cross-correlation formalism, contributing in tackling different astrophysical issues.– 24 – igure 9 . S/N from ∆ χ analysis (equation (5.2)) for discerning the two considered astrophysicalscenarios. Galaxies with ψ > , , M (cid:12) / yr are considered (top, central and bottom panel re-spectively). Continuous/dashed lines refer to the benchmark/Z-dependent model assumed as fiducial.Colors refer to different values of f sky as shown in legend. – 25 – Conclusions
In this work we have expanded the investigation in the field of the cross-correlations betweenresolved GWs signals and LSS tracers. We worked in the harmonic space with the numbercounts angular power spectra. The two categories of tracers we considered consist in resolvedGWs events from BH-BH, BH-NS, NS-NS mergers detectable by the Einstein Telescope andactively star-forming galaxies with SFR cuts of ψ > , , M (cid:12) / yr. We characterizedthem with their redshift distributions, bias and magnification bias values, presenting a detaileddescription of the computation of these quantities. We stress again that both the SFRF andthe GWs distributions derive from the same type of sources (galaxies) but trace them in adifferent way, since GW signals depend not only on the galaxy SFR, but also on the galaxymetallicity and on stellar and binary evolution prescriptions. For this reasons we kept intoaccount all the aforementioned elements in this work. Cross-correlating the same sources viatwo different messengers can help not only in alleviating systematics but also in enhancingthe amount of astrophysical information encoded in the signal.In our analysis we took into account all lensing and all general relativistic contributionsin the computation of the observed number counts fluctuations and we extended the basisfor future works regarding the GW × LSS cross-correlations with a more robust theoreticalastrophysical background.After computing the number counts angular power spectra for all the combinations of our GW × LSS tracers, we estimated Signal-to-Noise ratios in order to forecast the detectabilityof the cross-correlation signal. We have considered both tomographic and non-tomographicapproaches. Our results show that in several scenarios it is possible to reach a Signal-to-Noise ratio higher than unity, whereas it is not always the case for cross-correlations betweendistant redshift bins (in the tomographic case) or with a low number of observed objects.In addition, the total cumulative Signal-to-Noise ratios for each of the probes considered inthis work are in turn quite optimistic. Even though each of the many considerable surveys(such as ALMA [123], JWST [124], EMU [125], SKA [126] and many others) will have its ownspecifics, this work still provides a general possibility to gauge the cross-correlations efficacy.Finally, we have investigated the possibility of exploiting the GW × LSS formalism tocompare and test possible scenarios in the astrophysical modeling of GWs events. In par-ticular we considered a proof-of-concept case in which we made use of a ∆ χ statistics tocompare the cross-correlation signal obtained by modeling the COs merging efficiency with aspecific metallicity dependency with respect to a benchmark signal obtained neglecting thisdependency. We have showed that in principle, given enough individual objects observed (i.e.enough observation time, observed fraction of the sky, etc.) a metallicity dependency fea-ture could be discerned from the flat benchmark case. This is another step in the promisingmulti-tracers field and towards its astrophysical applications for future works to come. Acknowledgments
This work has been partially supported by PRIN MIUR 2017 prot. 20173ML3WW 002,‘Opening the ALMA window on the cosmic evolution of gas, stars and supermassive blackholes’. A.L. acknowledges the MIUR grant ‘Finanziamento annuale individuale attivitá basedi ricerca’ and the EU H2020-MSCA-ITN-2019 Project 860744 ‘BiD4BEST: Big Data appli-cations for Black hole Evolution STudies’. M.V. and G.S. are supported by INDARK PD51INFN grant. M.V. is also supported by ASI-INAF grant n.2017-14-H.0. We are thankful– 26 –o Nicola Bellomo, José Luis Bernal and Alvise Raccanelli for critical reading and helpfulsuggestions on an earlier version of this manuscript. L.B. acknowledges Martyna Chruslinskafor helpful discussions. We are thankful to the anonymous referee for thoughtful evaluationand helpful suggestions given to improve our manuscript.
A Relativistic number counts
We provide here the full expression for the relativistic number counts effects written in equa-tion (2.9): ∆ den (cid:96) ( k, z ) = b X δ ( k, τ z ) j (cid:96) , ∆ vel (cid:96) ( k, z ) = k H j (cid:48)(cid:48) (cid:96) V ( k, τ z ) + (cid:20) ( f evo X − H k j (cid:96) + (cid:18) H (cid:48) H + 2 − s X r ( z ) H + 5 s X − f evo X (cid:19) j (cid:48) (cid:96) (cid:21) V ( k, τ z ) , ∆ len (cid:96) ( k, z ) = (cid:96) ( (cid:96) + 1) 2 − s X (cid:90) r ( z )0 dr r ( z ) − rr ( z ) r [Φ( k, τ z ) + Ψ( k, τ z )] j (cid:96) ( kr ) , ∆ gr (cid:96) ( k, z ) = (cid:20)(cid:18) H (cid:48) H + 2 − s X r ( z ) H + 5 s X − f evo X + 1 (cid:19) Ψ( k, τ z ) + ( − s X ) Φ( k, τ z ) + H − Φ (cid:48) ( k, τ z ) (cid:21) j (cid:96) ++ (cid:90) r ( z )0 dr − s X r ( z ) [Φ( k, τ ) + Ψ( k, τ )] j (cid:96) ( kr ) , + (cid:90) r ( z )0 dr (cid:18) H (cid:48) H + 2 − s X r ( z ) H + 5 s X − f evo X (cid:19) r ( z ) (cid:2) Φ (cid:48) ( k, τ ) + Ψ (cid:48) ( k, τ ) (cid:3) j (cid:96) ( kr ) . (A.1)The meaning of the physical quantities written above is the following: b X is the bias param-eter, s X is the magnification bias parameter, f evo X is the evolution bias parameter, r is theconformal distance on the light cone, τ = τ − r is the conformal time, τ z = τ − r ( z ) , j (cid:96) , j (cid:48) (cid:96) = dj (cid:96) dy , j (cid:48)(cid:48) (cid:96) = d j (cid:96) dy are the Bessel functions and their derivatives (evaluated at y = kr ( z ) ifnot explicitly stated), H is the conformal Hubble parameter, the prime symbol (cid:48) stands forderivatives with respect to conformal time, δ is the density contrast in the comoving gauge, V is the peculiar velocity, Φ and Ψ are Bardeen potentials. B Angular power spectra for the tomographic case
We provide here the angular power spectra C (cid:96) ’s for all combinations of tracers in the tomo-graphic approach described in section 4.1. Note that each subplot reports cross-correlationsof redshift bins z i − z j , but each z i,j can actually be different between different tracers ac-cording to table 1. Some cases present power spectra with downward spikes: this is becausewe plotted the absolute values of the C (cid:96) ’s and the spikes simply correspond to change-of-signmultipoles (the reader interested in the non-surprising possibility of having negative C (cid:96) ’s canread e.g., ref. [127]). – 27 – − − − − z − z z − z BH-BH BH-NS NS-NS z − z z − z z − z z − z z − z z − z − − − − | C ‘ | z − z z − z z − z z − z z − z z − z z − z z − z − − − − z − z z − z z − z z − z z − z z − z z − z z − z − − − − z − z z − z ‘ z − z z − z z − z z − z z − z z − z galaxies: ψ > M (cid:12) / yr Figure 10 . Angular power spectra C (cid:96) ’s (absolute values) for the cross-correlations tomographic casesbetween galaxies with ψ > M (cid:12) / yr with all three types of GWs signals (BH-BH in red, BH-NS inblue, NS-NS in yellow). – 28 – − − − − z − z z − z BH-BH BH-NS NS-NS z − z z − z z − z z − z z − z z − z − − − − | C ‘ | z − z z − z z − z z − z z − z z − z z − z z − z − − − − z − z z − z z − z z − z z − z z − z z − z z − z − − − − z − z z − z ‘ z − z z − z z − z z − z z − z z − z galaxies: ψ > M (cid:12) / yr Figure 11 . Angular power spectra C (cid:96) ’s (absolute values) for the cross-correlations tomographic casesbetween galaxies with ψ > M (cid:12) / yr with all three types of GWs signals (BH-BH in red, BH-NS inblue, NS-NS in yellow). – 29 – − − − z − z z − z BH-BH BH-NS NS-NS z − z z − z z − z z − z − − − | C ‘ | z − z z − z z − z z − z z − z z − z − − − z − z z − z ‘ z − z z − z z − z z − z galaxies: ψ > M (cid:12) / yr Figure 12 . Angular power spectra C (cid:96) ’s (absolute values) for the cross-correlations tomographic casesbetween galaxies with ψ > M (cid:12) / yr with all three types of GWs signals (BH-BH in red, BH-NS inblue, NS-NS in yellow). References [1] The
LIGO Scientific Collaboration and Virgo Collaboration , B. P. Abbott et al.,“Observation of Gravitational Waves from a Binary Black Hole Merger”, Phys. Rev. Lett. (Feb, 2016) 061102, arXiv:1602.03837 .[2] The
LIGO Scientific Collaboration and Virgo Collaboration , B. P. Abbott et al.,“Properties of the Binary Black Hole Merger GW150914”, Phys. Rev. Lett. (Jun, 2016)241102, arXiv:1602.03840 .[3] S. Ho, C. Hirata, N. Padmanabhan, U. Seljak, and N. Bahcall, “Correlation of CMB with – 30 – arge-scale structure. I. Integrated Sachs-Wolfe tomography and cosmological implications”,Physical Review D no. 4, (2008) 043519.[4] C. M. Hirata, S. Ho, N. Padmanabhan, U. Seljak, and N. A. Bahcall, “Correlation of CMBwith large-scale structure. II. Weak lensing”, Physical Review D no. 4, (2008) 043520.[5] The Herschel ATLAS , F. Bianchini et al., “Cross-correlation between the CMB lensingpotential measured by Planck and high-z sub-mm galaxies detected by the Herschel-ATLASsurvey”, Astrophys. J. no. 1, (2015) 64, arXiv:1410.4502 [astro-ph.CO] .[6] F. Bianchini and A. Lapi, “Cross-correlation between cosmological and astrophysical datasets:the Planck and Herschel case”, IAU Symp. (2014) 202–205.[7] F. Bianchini et al., “Toward a tomographic analysis of the cross-correlation between PlanckCMB lensing and H-ATLAS galaxies”, Astrophys. J. no. 1, (2016) 24, arXiv:1511.05116[astro-ph.CO] .[8] S. Mukherjee, B. D. Wandelt, and J. Silk, “Multimessenger tests of gravity with weakly lensedgravitational waves”, Phys. Rev. D (May, 2020) 103509.[9] K. Fang, A. Banerjee, E. Charles, and Y. Omori, “A Cross-Correlation Study of High-energyNeutrinos and Tracers of Large-Scale Structure”, The Astrophysical Journal no. 2, (2020)112.[10] H. J. Martínez, M. E. Merchán, C. A. Valotto, and D. G. Lambas, “Quasar-galaxy andAGN-galaxy cross-correlations”, The Astrophysical Journal no. 2, (1999) 558.[11] B. Jain, R. Scranton, and R. K. Sheth, “QuasarâĂŤgalaxy and galaxyâĂŤgalaxycross-correlations: model predictions with realistic galaxies”, Monthly Notices of the RoyalAstronomical Society no. 1, (10, 2003) 62–70.[12] X. Yang, H. J. Mo, F. C. van den Bosch, S. M. Weinmann, C. Li, and Y. P. Jing, “Thecross-correlation between galaxies and groups: probing the galaxy distribution in and arounddark matter haloes”, Monthly Notices of the Royal Astronomical Society no. 2, (09, 2005)711–726.[13] K. Paech, N. Hamaus, B. Hoyle, M. Costanzi, T. Giannantonio, S. Hagstotz, G. Sauerwein,and J. Weller, “Cross-correlation of galaxies and galaxy clusters in the Sloan Digital SkySurvey and the importance of non-Poissonian shot noise”, Monthly Notices of the RoyalAstronomical Society no. 3, (06, 2017) 2566–2577.[14] A. Raccanelli, E. D. Kovetz, S. Bird, I. Cholis, and J. B. Muñoz, “Determining the progenitorsof merging black-hole binaries”, Phys. Rev. D (Jul, 2016) 023516, arXiv:1605.01405 .[15] G. Scelfo, N. Bellomo, A. Raccanelli, S. Matarrese, and L. Verde, “GW × LSS: chasing theprogenitors of merging binary black holes”, JCAP no. 09, (Sep, 2018) 039–039, arXiv:1809.03528v1 .[16] T. Namikawa, A. Nishizawa, and A. Taruya, “Anisotropies of gravitational-wave standardsirens as a new cosmological probe without redshift information”, Physical review letters no. 12, (2016) 121302.[17] D. Alonso, G. Cusin, P. G. Ferreira, and C. Pitrou, “Detecting the anisotropic astrophysicalgravitational wave background in the presence of shot noise through cross-correlations”, arXiv:2002.02888 .[18] D. Bertacca, A. Ricciardone, N. Bellomo, et al., “Projection effects on the observed angularspectrum of the astrophysical stochastic gravitational wave background”, Physical Review D no. 10, (May, 2020) 103513, arXiv:1909.11627 .[19] F. Calore, A. Cuoco, T. Regimbau, S. Sachdev, and P. D. Serpico, “Cross-correlating galaxycatalogs and gravitational waves: a tomographic approach”, Physical Review Research no. 2, (2020) 023314. – 31 –
20] S. Camera and A. Nishizawa, “Beyond Concordance Cosmology with Magnification ofGravitational-Wave Standard Sirens”, Phys. Rev. Lett. (Apr, 2013) 151103, arXiv:1303.5446 .[21] S. Mukherjee, B. D. Wandelt, and J. Silk, “Probing the theory of gravity with gravitationallensing of gravitational waves and galaxy surveys”, Monthly Notices of the RoyalAstronomical Society no. 2, (03, 2020) 1956–1970.[22] S. Mukherjee, B. D. Wandelt, S. M. Nissanke, and A. Silvestri, “Accurate and precisionCosmology with redshift unknown gravitational wave sources”, arXiv:2007.02943[astro-ph.CO] .[23] L. Boco, A. Lapi, S. Goswami, F. Perrotta, C. Baccigalupi, and L. Danese, “Merging Rates ofCompact Binaries in Galaxies: Perspectives for Gravitational Wave Detections”, TheAstrophysical Journal no. 2, (Aug., 2019) 157, arXiv:1907.06841 [astro-ph.GA] .[24] B. Sathyaprakash et al., “Scientific objectives of Einstein Telescope”, Classical and QuantumGravity no. 12, (2012) 124013, arXiv:1206.0331 .[25] A. F. Heavens and A. N. Taylor, “A spherical harmonic analysis of redshift space”, MonthlyNotices of the Royal Astronomical Society no. 2, (07, 1995) 483–497, arXiv:astro-ph/9409027 .[26] A. S. Szalay, T. Matsubara, and S. D. Landy, “Redshift-Space Distortions of the CorrelationFunction in Wide-Angle Galaxy Surveys”, The Astrophysical Journal no. 1, (May, 1998)L1–L4, arXiv:astro-ph/9712007 .[27] T. Matsubara, “The Correlation Function in Redshift Space: General Formula withWide-Angle Effects and Cosmological Distortions”, The Astrophysical Journal no. 1,(May, 2000) 1–23, arXiv:astro-ph/9908056 .[28] P. Pápai and I. Szapudi, “Non-perturbative effects of geometry in wide-angle redshiftdistortions”, Monthly Notices of the Royal Astronomical Society no. 1, (2008) 292–296, arXiv:0802.2940 .[29] A. Raccanelli, L. Samushia, and W. J. Percival, “Simulating redshift-space distortions forgalaxy pairs with wide angular separation”, Monthly Notices of the Royal AstronomicalSociety no. 4, (12, 2010) 1525–1533, arXiv:1006.1652 .[30] D. Bertacca, R. Maartens, A. Raccanelli, and C. Clarkson, “Beyond the plane-parallel andNewtonian approach: wide-angle redshift distortions and convergence in general relativity”,Journal of Cosmology and Astroparticle Physics no. 10, (Oct, 2012) 025, arXiv:1205.5221 .[31] A. Raccanelli, D. Bertacca, D. Pietrobon, F. Schmidt, L. Samushia, N. Bartolo, O. DorÃľ,S. Matarrese, and W. J. Percival, “Testing gravity using large-scale redshift-space distortions”,Monthly Notices of the Royal Astronomical Society no. 1, (09, 2013) 89–100, arXiv:1207.0500 .[32] L. Dai, M. Kamionkowski, and D. Jeong, “Total angular momentum waves for scalar, vector,and tensor fields”, Phys. Rev. D (Dec, 2012) 125013, arXiv:1209.0761 .[33] J. Yoo and V. Desjacques, “All-sky analysis of the general relativistic galaxy power spectrum”,Phys. Rev. D (Jul, 2013) 023502, arXiv:1301.4501 .[34] C. Blake, P. Carter, and J. Koda, “Power spectrum multipoles on the curved sky: anapplication to the 6-degree Field Galaxy Survey”, Monthly Notices of the Royal AstronomicalSociety no. 4, (2018) 5168–5183, arXiv:1801.04969 .[35] A. Taruya, S. Saga, M.-A. Breton, Y. Rasera, and T. Fujita, “Wide-angle redshift-spacedistortions at quasi-linear scales: cross-correlation functions from Zel’dovich approximation”, – 32 – onthly Notices of the Royal Astronomical Society no. 3, (11, 2019) 4162–4179, arXiv:1908.03854 .[36] A. Raccanelli, A. Bonaldi, M. Negrello, S. Matarrese, G. Tormen, and G. De Zotti, “Areassessment of the evidence of the Integrated Sachs-Wolfe effect through the WMAP-NVSScorrelation”, Monthly Notices of the Royal Astronomical Society no. 4, (2008) 2161–2166, arXiv:0802.0084 .[37] A. R. Pullen, T.-C. Chang, O. Doré, and A. Lidz, “Cross-correlations as a CosmologicalCarbon Monoxide Detector”, The Astrophysical Journal no. 1, (2013) 15, arXiv:1211.1397 .[38] C. Bonvin and R. Durrer, “What galaxy surveys really measure”, Phys. Rev. D (Sep, 2011)063505, arXiv:1105.5280 .[39] D. Blas, J. Lesgourgues, and T. Tram, “The Cosmic Linear Anisotropy Solving System(CLASS). Part II: Approximation schemes”, Journal of Cosmology and Astroparticle Physics no. 07, (2011) 034, arXiv:1104.2933 .[40] E. D. Dio, F. Montanari, J. Lesgourgues, and R. Durrer, “The CLASSgal code for relativisticcosmological large scale structure”, Journal of Cosmology and Astroparticle Physics no. 11, (2013) 044, arXiv:1307.1459 .[41] A. Challinor and A. Lewis, “Linear power spectrum of observed source number counts”, Phys.Rev. D (Aug, 2011) 043516, arXiv:1105.5292 .[42] N. Bellomo, J. L. Bernal, G. Scelfo, A. Raccanelli, and L. Verde, “Beware of commonly usedapproximations I: errors in forecasts”, arXiv:2005.10384 .[43] J. L. Bernal, N. Bellomo, A. Raccanelli, and L. Verde, “Beware of commonly usedapproximations II: estimating biases in the best-fit parameters”, arXiv:2005.09666 .[44] Planck Collaboration, Ade, P. A. R., et al., “Planck 2015 results - XIII. Cosmologicalparameters”, A&A (2016) A13, arXiv:1502.01589 .[45] N. Kaiser, “On the spatial correlations of Abell clusters”, The Astrophysical Journal (1984) L9–L12.[46] J. M. Bardeen, J. Bond, N. Kaiser, and A. Szalay, “The statistics of peaks of Gaussianrandom fields”, The Astrophysical Journal (1986) 15–61.[47] H. J. Mo and S. D. M. White, “An analytic model for the spatial clustering of dark matterhaloes”, Monthly Notices of the Royal Astronomical Society no. 2, (1996) 347–361, arXiv:astro-ph/9512127 .[48] S. Matarrese, P. Coles, F. Lucchin, and L. Moscardini, “Redshift evolution of clustering”,Monthly Notices of the Royal Astronomical Society no. 1, (03, 1997) 115–132, arXiv:astro-ph/9608004 .[49] A. Dekel and O. Lahav, “Stochastic Nonlinear Galaxy Biasing”, The Astrophysical Journal no. 1, (Jul, 1999) 24–34, arXiv:astro-ph/9806193 .[50] A. J. Benson, S. Cole, C. S. Frenk, C. M. Baugh, and C. G. Lacey, “The nature of galaxy biasand clustering”, Monthly Notices of the Royal Astronomical Society no. 4, (02, 2000)793–808, arXiv:astro-ph/9903343 .[51] J. A. Peacock and R. E. Smith, “Halo occupation numbers and galaxy bias”, Monthly Noticesof the Royal Astronomical Society no. 4, (11, 2000) 1144–1156, arXiv:astro-ph/0005010 .[52] V. Desjacques, D. Jeong, and F. Schmidt, “Large-scale galaxy bias”, Physics Reports (2018) 1 – 193. – 33 –
53] E. L. Turner, J. P. Ostriker, and J. R. Gott, III, “The statistics of gravitational lenses - Thedistributions of image angular separations and lens redshifts”, Astrophysical Journal (Sep, 1984) 1–22.[54] A. Challinor and A. Lewis, “Linear power spectrum of observed source number counts”, Phys.Rev. D (Aug, 2011) 043516, arXiv:1105.5292 .[55] D. Jeong, F. Schmidt, and C. M. Hirata, “Large-scale clustering of galaxies in generalrelativity”, Phys. Rev. D (Jan, 2012) 023504, arXiv:1107.5427 .[56] D. Bertacca, R. Maartens, A. Raccanelli, and C. Clarkson, “Beyond the plane-parallel andNewtonian approach: wide-angle redshift distortions and convergence in general relativity”,Journal of Cosmology and Astroparticle Physics no. 10, (2012) 025, arXiv:1205.5221 .[57] G. R. Meurer, T. M. Heckman, and D. Calzetti, “Dust Absorption and the UltravioletLuminosity Density at z ~3 as Calibrated by Local Starburst Galaxies”, The AstrophysicalJournal no. 1, (Aug., 1999) 64–80, arXiv:astro-ph/9903054 [astro-ph] .[58] D. Calzetti, L. Armus, R. C. Bohlin, A. L. Kinney, J. Koornneef, and T. Storchi-Bergmann,“The Dust Content and Opacity of Actively Star-forming Galaxies”, The AstrophysicalJournal no. 2, (Apr., 2000) 682–695, arXiv:astro-ph/9911459 [astro-ph] .[59] R. J. Bouwens, G. D. Illingworth, P. A. Oesch, et al., “UV Luminosity Functions at Redshiftsz ∼ ∼
10: 10,000 Galaxies from HST Legacy Fields”, The Astrophysical Journal no. 1, (Apr., 2015) 34, arXiv:1403.4295 [astro-ph.CO] .[60] L. Silva, G. L. Granato, A. Bressan, and L. Danese, “Modeling the Effects of Dust on GalacticSpectral Energy Distributions from the Ultraviolet to the Millimeter Band”, TheAstrophysical Journal no. 1, (Dec., 1998) 103–117.[61] A. Efstathiou, M. Rowan-Robinson, and R. Siebenmorgen, “Massive star formation ingalaxies: radiative transfer models of the UV to millimetre emission of starburst galaxies”,Monthly Notices of the Royal Astronomical Society no. 4, (Apr., 2000) 734–744, arXiv:astro-ph/9912252 [astro-ph] .[62] K. E. K. Coppin, J. E. Geach, O. Almaini, et al., “The SCUBA-2 Cosmology Legacy Survey:the submillimetre properties of Lyman-break galaxies at z = 3-5”, Monthly Notices of theRoyal Astronomical Society no. 2, (Jan., 2015) 1293–1304, arXiv:1407.6712[astro-ph.GA] .[63] N. A. Reddy, M. Kriek, A. E. Shapley, et al., “The MOSDEF Survey: Measurements ofBalmer Decrements and the Dust Attenuation Curve at Redshifts z ~1.4-2.6”, TheAstrophysical Journal no. 2, (June, 2015) 259, arXiv:1504.02782 [astro-ph.GA] .[64] Y. Fudamoto, P. A. Oesch, E. Schinnerer, et al., “The dust attenuation of star-forminggalaxies at z ∼ no. 1, (Nov., 2017) 483–490, arXiv:1705.01559[astro-ph.GA] .[65] A. Lapi, J. González-Nuevo, L. Fan, et al., “Herschel-ATLAS Galaxy Counts andHigh-redshift Luminosity Functions: The Formation of Massive Early-type Galaxies”, TheAstrophysical Journal no. 1, (Nov., 2011) 24, arXiv:1108.3911 [astro-ph.CO] .[66] C. Gruppioni, F. Pozzi, G. Rodighiero, et al., “The Herschel PEP/HerMES luminosityfunction - I. Probing the evolution of PACS selected Galaxies to z (cid:39) no. 1, (June, 2013) 23–52, arXiv:1302.5209[astro-ph.CO] .[67] C. Gruppioni, F. Calura, F. Pozzi, et al., “Star formation in Herschel’s Monsters versussemi-analytic models”, Monthly Notices of the Royal Astronomical Society no. 4, (Aug.,2015) 3419–3426, arXiv:1506.01518 [astro-ph.GA] . – 34 –
68] C. Gruppioni and F. Pozzi, “On the existence of bright IR galaxies at z > 2: tensionbetween Herschel and SCUBA-2 results?”, Monthly Notices of the Royal Astronomical Society no. 2, (Feb., 2019) 1993–1999, arXiv:1812.00682 [astro-ph.GA] .[69] B. Magnelli, P. Popesso, S. Berta, et al., “The deepest Herschel-PACS far-infrared survey:number counts and infrared luminosity functions from combined PEP/GOODS-Hobservations”, A&A (May, 2013) A132, arXiv:1303.4436 [astro-ph.CO] .[70] M. Rowan-Robinson, S. Oliver, L. Wang, D. Farrah, D. L. Clements, C. Gruppioni,L. Marchetti, D. Rigopoulou, and M. Vaccari, “The star formation rate density from z = 1 to6”, Monthly Notices of the Royal Astronomical Society no. 1, (Sept., 2016) 1100–1111, arXiv:1605.03937 [astro-ph.GA] .[71] J. S. Dunlop, R. J. McLure, A. D. Biggs, et al., “A deep ALMA image of the Hubble UltraDeep Field”, Monthly Notices of the Royal Astronomical Society no. 1, (Apr., 2017)861–883, arXiv:1606.00227 [astro-ph.GA] .[72] D. Liu, E. Daddi, M. Dickinson, et al., ““Super-deblended” Dust Emission in Galaxies. I. TheGOODS-North Catalog and the Cosmic Star Formation Rate Density out to Redshift 6”, TheAstrophysical Journal no. 2, (Feb., 2018) 172, arXiv:1703.05281 [astro-ph.GA] .[73] D. A. Riechers, T. K. D. Leung, R. J. Ivison, et al., “Rise of the Titans: A Dusty,Hyper-luminous “870 µ m Riser” Galaxy at z ∼ no. 1,(Nov., 2017) 1, arXiv:1705.09660 [astro-ph.GA] .[74] D. P. Marrone, J. S. Spilker, C. C. Hayward, et al., “Galaxy growth in a massive halo in thefirst billion years of cosmic history”, Nature no. 7686, (Jan., 2018) 51–54, arXiv:1712.03020 [astro-ph.GA] .[75] J. A. Zavala, A. Montaña, D. H. Hughes, et al., “A dusty star-forming galaxy at z = 6revealed by strong gravitational lensing”, Nature Astronomy (Nov., 2018) 56–62, arXiv:1707.09022 [astro-ph.GA] .[76] M. Novak, V. Smolčić, J. Delhaize, et al., “The VLA-COSMOS 3 GHz Large Project: Cosmicstar formation history since z ∼ (2017) A5.[77] C. Mancuso, A. Lapi, J. Shi, J. Gonzalez-Nuevo, R. Aversa, and L. Danese, “The Quest forDusty Star-forming Galaxies at High Redshift z (cid:38) no. 2,(June, 2016) 128, arXiv:1604.02507 [astro-ph.GA] .[78] R. Aversa, A. Lapi, G. de Zotti, F. Shankar, and L. Danese, “Black Hole and GalaxyCoevolution from Continuity Equation and Abundance Matching”, The Astrophysical Journal no. 1, (Sept., 2015) 74, arXiv:1507.07318 [astro-ph.GA] .[79] A. Lapi, F. Shankar, J. Mao, G. Granato, L. Silva, G. De Zotti, and L. Danese, “Quasarluminosity functions from joint evolution of black holes and host galaxies”, The AstrophysicalJournal no. 1, (2006) 42.[80] A. Lapi and A. Cavaliere, “Self-similar Dynamical Relaxation of Dark Matter Halos in anExpanding Universe”, The Astrophysical Journal no. 2, (2011) 127.[81] A. Lapi, S. Raimundo, R. Aversa, Z.-Y. Cai, M. Negrello, A. Celotti, G. De Zotti, andL. Danese, “The coevolution of supermassive black holes and massive galaxies at highredshift”, The Astrophysical Journal no. 2, (2014) 69.[82] R. K. Sheth, H. J. Mo, and G. Tormen, “Ellipsoidal collapse and an improved model for thenumber and spatial distribution of dark matter haloes”, Monthly Notices of the RoyalAstronomical Society no. 1, (May, 2001) 1–12, arXiv:astro-ph/9907024 [astro-ph] .[83] A. Lapi and L. Danese, “Statistics of dark matter halos in the excursion set peak framework”,JCAP no. 7, (July, 2014) 044, arXiv:1407.1137 [astro-ph.CO] . – 35 –
84] F. Bernardeau, S. Colombi, E. GaztaÃśaga, and R. Scoccimarro, “Large-scale structure of theUniverse and cosmological perturbation theory”, Physics Reports no. 1, (2002) 1 – 248.[85] M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mand el, andR. O’Shaughnessy, “Double Compact Objects. II. Cosmological Merger Rates”, TheAstrophysical Journal no. 1, (Dec., 2013) 72, arXiv:1308.1546 [astro-ph.HE] .[86] M. Dominik, E. Berti, R. O’Shaughnessy, I. Mandel, K. Belczynski, C. Fryer, D. E. Holz,T. Bulik, and F. Pannarale, “Double Compact Objects III: Gravitational-wave DetectionRates”, The Astrophysical Journal no. 2, (June, 2015) 263, arXiv:1405.7016[astro-ph.HE] .[87] S. E. de Mink, N. Langer, R. G. Izzard, H. Sana, and A. de Koter, “The Rotation Rates ofMassive Stars: The Role of Binary Interaction through Tides, Mass Transfer, and Mergers”,The Astrophysical Journal no. 2, (Feb., 2013) 166, arXiv:1211.3742 [astro-ph.SR] .[88] M. Spera, M. Mapelli, and A. Bressan, “The mass spectrum of compact remnants from thePARSEC stellar evolution tracks”, Monthly Notices of the Royal Astronomical Society no. 4, (Aug., 2015) 4086–4103, arXiv:1505.05201 [astro-ph.SR] .[89] M. Spera and M. Mapelli, “Very massive stars, pair-instability supernovae andintermediate-mass black holes with the sevn code”, Monthly Notices of the RoyalAstronomical Society no. 4, (Oct., 2017) 4739–4749, arXiv:1706.06109 [astro-ph.SR] .[90] M. Spera, M. Mapelli, N. Giacobbo, A. A. Trani, A. Bressan, and G. Costa, “Merging blackhole binaries with the SEVN code”, Monthly Notices of the Royal Astronomical Society no. 1, (May, 2019) 889–907, arXiv:1809.04605 [astro-ph.HE] .[91] N. Giacobbo and M. Mapelli, “The progenitors of compact-object binaries: impact ofmetallicity, common envelope and natal kicks”, Monthly Notices of the Royal AstronomicalSociety no. 2, (Oct., 2018) 2011–2030, arXiv:1806.00001 [astro-ph.HE] .[92] M. Mapelli, N. Giacobbo, E. Ripamonti, and M. Spera, “The cosmic merger rate of stellarblack hole binaries from the Illustris simulation”, Monthly Notices of the Royal AstronomicalSociety no. 2, (Dec., 2017) 2422–2435, arXiv:1708.05722 [astro-ph.GA] .[93] R. O’Shaughnessy, J. M. Bellovary, A. Brooks, S. Shen, F. Governato, and C. R. Christensen,“The effects of host galaxy properties on merging compact binaries detectable by LIGO”,Monthly Notices of the Royal Astronomical Society no. 3, (Jan., 2017) 2831–2839, arXiv:1609.06715 [astro-ph.GA] .[94] A. Lamberts, S. Garrison-Kimmel, P. F. Hopkins, E. Quataert, J. S. Bullock, C. A.Faucher-Giguère, A. Wetzel, D. Kereš, K. Drango, and R. E. Sand erson, “Predicting thebinary black hole population of the Milky Way with cosmological simulations”, MonthlyNotices of the Royal Astronomical Society no. 2, (Oct., 2018) 2704–2718, arXiv:1801.03099 [astro-ph.GA] .[95] M. Mapelli and N. Giacobbo, “The cosmic merger rate of neutron stars and black holes”,Monthly Notices of the Royal Astronomical Society no. 4, (Oct., 2018) 4391–4398, arXiv:1806.04866 [astro-ph.HE] .[96] M. C. Artale, M. Mapelli, N. Giacobbo, N. B. Sabha, M. Spera, F. Santoliquido, andA. Bressan, “Host galaxies of merging compact objects: mass, star formation rate, metallicity,and colours”, Monthly Notices of the Royal Astronomical Society no. 2, (Aug., 2019)1675–1688, arXiv:1903.00083 [astro-ph.GA] .[97] K. Belczynski, M. Dominik, T. Bulik, R. O’Shaughnessy, C. Fryer, and D. E. Holz, “TheEffect of Metallicity on the Detection Prospects for Gravitational Waves”, The AstrophysicalJournal Letters no. 2, (2010) L138, arXiv:1004.0386 .[98] A. Lamberts, S. Garrison-Kimmel, D. R. Clausen, and P. F. Hopkins, “When and where did – 36 –
W150914 form?”, Monthly Notices of the Royal Astronomical Society no. 1, (Nov.,2016) L31–L35, arXiv:1605.08783 [astro-ph.HE] .[99] L. Cao, Y. Lu, and Y. Zhao, “Host galaxy properties of mergers of stellar binary black holesand their implications for advanced LIGO gravitational wave sources”, Monthly Notices of theRoyal Astronomical Society no. 4, (Mar., 2018) 4997–5007, arXiv:1711.09190[astro-ph.GA] .[100] O. D. Elbert, J. S. Bullock, and M. Kaplinghat, “Counting black holes: The cosmic stellarremnant population and implications for LIGO”, Monthly Notices of the Royal AstronomicalSociety no. 1, (Jan., 2018) 1186–1194, arXiv:1703.02551 [astro-ph.GA] .[101] S.-S. Li, S. Mao, Y. Zhao, and Y. Lu, “Gravitational lensing of gravitational waves: astatistical perspective”, Monthly Notices of the Royal Astronomical Society no. 2, (May,2018) 2220–2229, arXiv:1802.05089 [astro-ph.CO] .[102] C. J. Neijssel, A. Vigna-Gómez, S. Stevenson, J. W. Barrett, S. M. Gaebel, F. S.Broekgaarden, S. E. de Mink, D. Szécsi, S. Vinciguerra, and I. Mandel, “The effect of themetallicity-specific star formation history on double compact object mergers”, MonthlyNotices of the Royal Astronomical Society no. 3, (Dec., 2019) 3740–3759, arXiv:1906.08136 [astro-ph.SR] .[103] L. Pantoni, A. Lapi, M. Massardi, S. Goswami, and L. Danese, “New Analytic Solutions forGalaxy Evolution: Gas, Stars, Metals, and Dust in Local ETGs and Their High-zStar-forming Progenitors”, The Astrophysical Journal no. 2, (Aug., 2019) 129, arXiv:1906.07458 [astro-ph.GA] .[104] A. Lapi, L. Pantoni, L. Boco, and L. Danese, “New Analytic Solutions for Galaxy EvolutionII: Wind Recycling, Galactic Fountains and Late-Type Galaxies”, arXiv e-prints (May, 2020)arXiv:2006.01643, arXiv:2006.01643 [astro-ph.GA] .[105] M. Arrigoni, S. C. Trager, R. S. Somerville, and B. K. Gibson, “Galactic chemical evolution inhierarchical formation models - I. Early-type galaxies in the local Universe”, Monthly Noticesof the Royal Astronomical Society no. 1, (Feb., 2010) 173–190, arXiv:0905.4189[astro-ph.CO] .[106] M. Spolaor, C. Kobayashi, D. A. Forbes, W. J. Couch, and G. K. T. Hau, “Early-typegalaxies at large galactocentric radii - II. Metallicity gradients and the [Z/H]-mass,[ α /Fe]-mass relations”, Monthly Notices of the Royal Astronomical Society no. 1, (Oct.,2010) 272–292, arXiv:1006.1698 [astro-ph.CO] .[107] A. Gallazzi, E. F. Bell, S. Zibetti, J. Brinchmann, and D. D. Kelson, “Charting the Evolutionof the Ages and Metallicities of Massive Galaxies since z = 0.7”, The Astrophysical Journal no. 1, (June, 2014) 72, arXiv:1404.5624 [astro-ph.GA] .[108] B. H. Andrews and P. Martini, “The Mass-Metallicity Relation with the Direct Method onStacked Spectra of SDSS Galaxies”, The Astrophysical Journal no. 2, (Mar., 2013) 140, arXiv:1211.3418 [astro-ph.CO] .[109] H. J. Zahid, D. Kashino, J. D. Silverman, et al., “The FMOS-COSMOS Survey ofStar-forming Galaxies at z ~1.6. II. The Mass-Metallicity Relation and the Dependence onStar Formation Rate and Dust Extinction”, The Astrophysical Journal no. 1, (Sept.,2014) 75, arXiv:1310.4950 [astro-ph.CO] .[110] I. G. de la Rosa, F. La Barbera, I. Ferreras, J. Sánchez Almeida, C. Dalla Vecchia,I. Martínez-Valpuesta, and M. Stringer, “The fate of high-redshift massive compact galaxies”,Monthly Notices of the Royal Astronomical Society no. 2, (Apr., 2016) 1916–1930, arXiv:1601.03920 [astro-ph.GA] .[111] M. Onodera, C. M. Carollo, S. Lilly, A. Renzini, et al., “ISM Excitation and Metallicity of – 37 – tar-forming Galaxies at z (cid:39) no. 1, (May, 2016) 42, arXiv:1602.02779 [astro-ph.GA] .[112] G. Chabrier, “The Galactic Disk Mass Function: Reconciliation of the Hubble SpaceTelescope and Nearby Determinations”, The Astrophysical Journal Letters no. 2, (Apr.,2003) L133–L136, arXiv:astro-ph/0302511 [astro-ph] .[113] E. D. Kovetz, I. Cholis, P. C. Breysse, and M. Kamionkowski, “Black hole mass function fromgravitational wave measurements”, Physical Review D no. 10, (2017) 103010.[114] M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mand el, andR. O’Shaughnessy, “Double Compact Objects. I. The Significance of the Common Envelopeon Merger Rates”, The Astrophysical Journal no. 1, (Nov., 2012) 52, arXiv:1202.4901[astro-ph.HE] .[115] S. E. de Mink and K. Belczynski, “Merger Rates of Double Neutron Stars and Stellar OriginBlack Holes: The Impact of Initial Conditions on Binary Evolution Predictions”, TheAstrophysical Journal no. 1, (Nov., 2015) 58, arXiv:1506.03573 [astro-ph.HE] .[116] M. Chruslinska, G. Nelemans, and K. Belczynski, “The influence of the distribution of cosmicstar formation at different metallicities on the properties of merging double compact objects”,Monthly Notices of the Royal Astronomical Society no. 4, (Feb., 2019) 5012–5017, arXiv:1811.03565 [astro-ph.HE] .[117] R. O’Shaughnessy, V. Kalogera, and K. Belczynski, “Binary Compact Object CoalescenceRates: The Role of Elliptical Galaxies”, The Astrophysical Journal no. 1, (June, 2010)615–633, arXiv:0908.3635 [astro-ph.CO] .[118] M. Chruslinska, K. Belczynski, J. Klencki, and M. Benacquista, “Double neutron stars:merger rates revisited”, Monthly Notices of the Royal Astronomical Society no. 3, (Mar.,2018) 2937–2958, arXiv:1708.07885 [astro-ph.HE] .[119] The LIGO Scientific Collaboration and Virgo Collaboration , Abbott, B. P. andothers, “GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary MergersObserved by LIGO and Virgo during the First and Second Observing Runs”, Physical ReviewX no. 3, (July, 2019) 031040, arXiv:1811.12907 [astro-ph.HE] .[120] S. R. Taylor and J. R. Gair, “Cosmology with the lights off: Standard sirens in the EinsteinTelescope era”, Physical Review D no. 2, (July, 2012) 023502, arXiv:1204.6739[astro-ph.CO] .[121] L. Boco, A. Lapi, and L. Danese, “Growth of Supermassive Black Hole Seeds in ETGStar-forming Progenitors: Multiple Merging of Stellar Compact Remnants via GaseousDynamical Friction and Gravitational-wave Emission”, The Astrophysical Journal no. 1,(Mar., 2020) 94, arXiv:2002.03645 [astro-ph.GA] .[122] S. Klimenko, G. Vedovato, Drago, and others., “Localization of gravitational wave sourceswith networks of advanced detectors”, Physical Review D no. 10, (May, 2011) 102001, arXiv:1101.5408 [astro-ph.IM] .[123] A. Wootten and A. R. Thompson, “The Atacama Large Millimeter/Submillimeter Array”,Proceedings of the IEEE no. 8, (2009) 1463–1471.[124] J. P. Gardner et al., “The james webb space telescope”, Space Science Reviews no. 4,(2006) 485–606.[125] R. P. Norris et al., “EMU: Evolutionary Map of the Universe”, Publications of theAstronomical Society of Australia no. 3, (2011) 215–248, arXiv:1106.3219 .[126] R. Maartens et al., “Cosmology with the SKA - overview”, arXiv:1501.04076 . – 38 – no. 01, (2014) 042.no. 01, (2014) 042.