Magnification and evolution bias of transient sources: GWs and SNIa

Third-generation gravitational wave (GW) observatories such as the Einstein Telescope and Cosmic Explorer, together with the LSST survey at the Vera Rubin Observatory, will yield an abundance of extra-galactic transient objects. This opens the exciting possibility of using GW sources and Supernovae Type Ia (SNIa) as luminosity distance tracers of large-scale structure for the first time. The large volumes accessible to these surveys imply that we may need to include relativistic corrections, such as lensing and Doppler magnification. However, the amplitude of these effects depends on the magnification and evolution biases of the transient sources, which are not yet understood. In this paper we develop comprehensive frameworks to address and model these biases for both populations of transient objects; in particular, we define how to compute these biases for GW sources. We then analyse the impact of magnification and evolution biases on the relativistic corrections and on the angular power spectrum of these sources. We show that correct modelling and implementation of these biases is crucial for measuring the cross-correlations of transient sources at higher redshifts.

Along with the sheer number of objects seen, this new generation of observatories will push the horizon of detections further: whilst LSST will observe SNIa events at z < 4 [21,22,28], ET and CE are predicted to access virtually all binary black hole (BBH) mergers up to z ∼ 10 [28, [30][31][32]35].Greatly increasing the source distances also increases the relevance of relativistic correction effects.This is due to both the impact of lensing on the larger redshifts covered, and due to tracers evolving with cosmic time, thus not distributing across different redshift bins equally [18,36].
These future prospects indicate that the kinds of clustering analyses carried out with galaxies and IM will soon be applicable to transient events such as SNIa and GWs too.A key subtlety is that the former are tracers living in redshift space, whilst the latter carry distance information only under the form of a luminosity distance.In a previous work [37] we explored the difference between clustering analysis in redshift space and in luminosity distance space for a generic tracer.We found that the two can be significantly different, leading to large discrepancies in the corresponding angular power spectrum, up to 50% at large scales.Therefore, any analysis utilising tracers such as GWs or SNIa should be carried out in luminosity distance.
Primarily, one is required to build an expression for the observed number counts in luminosity distance space, which not only traces the underlying density of matter, but also includes distortion effects along our past light-cone.One can write the generic expression for the density contrast in luminosity distance space as [37,38]: where γ ≡ rH/1 + rH with comoving distance r and comoving Hubble parameter H, δ n is the density contrast in the Newtonian gauge, D L the luminosity distance, and V the volume.
The distortion is then found to be dependent on two extra functions, the evolution bias [39] b e ≡ ∂ ln n ∂ ln a d th , where d th is the detector's observation threshold (taken as the luminosity cut L c in the standard galaxy case and signal-to-noise ratio threshold ρ th for GWs), and the magnification bias s, defined as the change in the comoving number density n at fixed redshift/distance with respect to the luminosity cut.Instead, b e is the change in the comoving number density with respect to the scale factor, while keeping the detector's threshold fixed [39].Physically, the magnification bias accounts for objects being magnified in or out of the detector's flux limit as a result of a perturbation, and the evolution bias describes the impact on clustering of a (possibly) non-conserved comoving number density through cosmic time.A null evolution bias would imply a constant observed number density through redshift, whilst a non-zero one corrects each redshift bin for the effect of an evolving population.In other words, it traces how the observed number of objects per unit volume changes as the universe expands.We define the magnification bias for GWs and SNIa differently in order to preserve the general expression in eq.(1.1): the former in terms of signal-to-noise ratio (SNR) threshold and the latter in terms of magnitude cut, respectively.In the following sections we will provide the expression for each tracer, as well as specify the magnification bias for different surveys and thresholds.
Here, we propose a theoretical framework to model both GWs and SNIa, thus aiming at consistently translating terminology and modelling previously used in galaxy clustering into equivalent expressions applicable to these new transient tracers.One should note that some of these biases, at least for GWs, were recently explored for the first time in [24][25][26]40], particularly in the context of cross-correlations, as these have become increasingly exciting with the prospect of future detectors.
The paper is structured as follows.In section 2 we describe the modelling of the biases for GWs, justifying our theoretical formalism and producing the required event rate and chirp mass distributions.Further, we model the biases for SNIa in section 3. Finally, we apply these biases to the kernels of the number count fluctuation and the angular power spectrum in section 4, and study their impact.Section 5 is then devoted to summary and conclusions.

GW Biases
Here we illustrate the modelling of the magnification and evolution biases for GWs.We first explain the reasoning behind choosing the signal-to-noise ratio instead of luminosity as discriminant for the magnification bias; we then describe both the modelling for the event rate and the chirp mass distribution, which we derive using the primary mass distribution from [41][42][43].An interested reader can find the distributions of primary and secondary masses in Appendix A, and the derivation of the chirp mass probability distribution function in Appendix B. We then evaluate the GW biases for the different mass models used.
For galaxy surveys the magnification bias is generally defined from [39] as: where L c is the luminosity threshold of our survey at each redshift and n g (a, L c ) is the comoving number density of sources above the threshold.The latter is defined by integrating the (comoving) luminosity function Φ over luminosity: where L * is a characteristic luminosity in the luminosity function.The number of sources that are above the corresponding luminosity cut is the same as the number of sources that are observed.This implies that regardless of the choice of luminosity cut, n g can be zero at a certain redshift, as there will be no sources emitting at a sufficiently high luminosity.Therefore, the value for the magnification bias is expected to increase with redshift up until its validity limit, i.e.where sources are not detectable anymore, meaning that observing fewer and fewer sources leads to a (positive) larger value of s.
We can now transport this concept from the galaxy case to the GW scenario.However, for GWs we can't produce a detector-independent quantity analogous to EM luminosity, as the quality of a detection (the SNR, ρ) is inherently tied to the one-sided power spectral density of the detector, i.e. its sensitivity, P SD(f ) [44]: In galaxy surveys the samples are also telescope dependent, as different telescopes target different parts of the spectrum or optimize for different types of galaxies.As an example Euclid will target Hα emitting galaxies, while DESI will have broader target to any emission line.The response function of a GW detector is equivalent to the instrumental design of optical and near-infrared telescopes.Using the characteristic strain h(f ) of a GW event (which depends on chirp mass M and redshift) is not possible either, as it still has to be compared to the detector's sensitivity curve, and their ratio integrated to determine whether the event is detected and at which level of confidence.The two scenarios would be similar if we assumed the sensitivity curve of GWs detector to be flat across all frequencies, and the GWs signal to be a single-frequency burst.All signals would then see the same response from the detector, as they would not move across frequency space.We therefore define the magnification bias for GW merger events in term of the SNR as The factor of 1 5 comes from implicitly fixing the expression eq.(1.1) and setting s accordingly; this was done to keep eq.(1.1) general, particularly for coding purposes.We therefore need to model the number density of observed events as a function of the detector signal-to-noise ratio ρ.Eq. (2.4) defines the magnification bias for any GW source.However, for the purpose of this paper we will only focus on GWs from binary black hole mergers.The same method can be applied to binary neutron stars and neutron star black hole pairs, if the appropriate number density is provided.

Event rate for GWs
The number density of observed sources is the number above a certain SNR threshold, which is usually assumed to be ρ th = 8 for a single detector.In the case of multiple detectors, ρ th is added in quadrature.
From [45], we can model the number density of BBH mergers in comoving volume as: where τ is the observation time of the detector, R GW (z) is the intrinsic merger rate, ϕ(M) is the chirp mass M distribution, with Selection effects of the signal-to-noise threshold are included through the implementation of a survival function S(ρ th ; M, z), which is defined as the fraction of sources that are above the SNR threshold, ρ th , at any given mass and redshift bin1 .We note that for physical experiments, the sensitivity curve should be defined on the detected SNR, rather than the optimal one.The former has a non-central χ 2 distribution [46] which depends on the optimal SNR, which by itself is a function of the PSD.Although our simplified assumption is not expected to qualitatively change the result of this study, a rigorous analysis with real data would have to take this into account for a more realistic experiment.To evaluate S(ρ th ; M, z), we start from an expression for the SNR ρ for a compact binary merger [44,45]: where θ is the orientation of the binary with respect to the detector, and ρ 0 (M, z) encapsulates signal and detector's response.For simplicity, here we only use a 0PN approximation for the signal to compute the SNR.A more rigorous analysis would include higher order corrections to describe merger and ringdown parts without this approximation, thus including a boost to the SNR.Assuming random orientations, the PDF of θ can be well approximated by [44]: Left: Fraction of sources detected, S(ρ th /ρ 0 ) as a function of ρ 0 , i.e. the characteristic SNR of the source.Low values of S imply either a low SNR threshold or loud events, whereas large S signifies a high threshold or quiet events.Center: Fraction of sources undetected at each random orientation.θ is equivalent to the ratio between the SNR threshold ρ th and the characteristic SNR of the source, ρ 0 .Right: Orientation PDF.
for 0 < θ < 4 and P (θ) = 0 otherwise.Note that this is valid for a single L-shaped detector.In the case of three detectors such as LVK, the SNR can be added in quadrature, although for a triangular configuration such as ET, this would change.However, for simplicity we assume the same function applies to an ET-like experiment.Further, we assume that f max corresponds to the frequency at the innermost stable circular orbit: From eq. (2.7) we have θ = ρ/ρ 0 , thus by fixing ρ = ρ th we can obtain a PDF for the sources that produce a sufficiently high SNR to be detected.The corresponding cumulative distribution function is the integral of eq.(2.9), evaluated from θ c to 4. Thus, the survival function S(ρ th ; M, z) from eq. (2.5) is defined as: where T (θ) is the integral of the angular orientation PDF P (θ). Figure 1 shows plots for these functions.Combining these results, the expression for the number density of sources becomes assuming an observation time τ and a redshift bin z.While the survival function is detector dependent, the chirp mass distribution depends on the mass distribution of black holes.To fully compute the biases we need to specify how likely are chirp masses of value M.
The full derivation of the chirp mass PDF can be found in Appendix B2 .Here we report only the final result for the PDF ϕ(M): where f and g are the distributions of secondary and primary masses respectively.x 1 is a real solution for the secondary mass m 2 in terms of chirp M and primary m 1 masses, given by the expression Note that this expression is obtained by solving the cubic equation (B.1) and is valid for We then employ current phenomenological distributions of primary and secondary masses from the catalogues of the LIGO-Virgo-KAGRA (LVK) collaboration [42,43], to describe f and g in eq.(2.13).Their full prescriptions can be found in Appendix A.

Magnification and evolution biases for GW events
Using the definition of magnification bias described in eq. ( 2.1), we can write: Initially, we plug in the sensitivity curve for a network of aLIGO-like interferometers to mimic the LVK detectors3 , adding the SNR in quadrature for three detectors.We find that the values of the magnification bias in this case are strongly dependent on the low number of sources detectable.In fact, figure 2 shows very large values of s for all models considered.Clearly, assuming a higher threshold of detection drastically decreases the number of observed sources even at lower redshift, and, consequently, increases the magnification bias.Conversely, a lower value of ρ th (e.g.blue line) yields lower values of s across a larger redshift range.Further, it is interesting to note that, at the same detector's threshold, different distributions of chirp masses have similar magnification biases.This suggests that for the chirp mass models assumed in this paper, s does not vary significantly.We note that this conclusion is not necessarily true for any chirp mass distribution.
The large values of s for LVK shown in figure 2 can be further explained with figure 3, investigating the dependence of the observed number density of GWs as seen by a LVKlike experiment on both redshift and SNR threshold.We calculate the observed number densities using (2.12), assuming an observation time of 1 year and that the intrinsic merger rate R GW (z) follows directly the Madau-Dickinson rate [47,48]:  Magnification bias for an LVK-era experiment.We stop calculating s for each model when the number of observed sources is below the arbitrary value of 100 events Gpc −3 ; this is to suggest a limit of validity for the bias, as clustering analysis requires a larger number of sources.Solid lines indicate a Power Law + Peak model for the primary BH mass, whilst dashed ones show the Broken Power Law model.Fainter lines represent the range where the biases become more uncertain, as the number of observed sources crosses 120, getting closer to our arbitrary cut.
with R 0 providing the merger rate at z = 0 and is given by [42,43] as R 0 = 23.9Gpc −3 yr −1 .In particular, the left plot shows how a larger value of ρ th will yield a steeper curve at low redshift, and, conversely, a lower value of ρ th produces a more gentle one.Therefore, as redshift increases, so does the distance between each curve, showing the number of sources missed/detected when changing the minimum SNR.This is made more explicit in the plot on the right, where curves at fixed redshifts are plotted against ρ th .Notably, the magnification bias can be read off directly here, as it is the slope at a given value of SNR threshold for a fixed redshift.Thus, one can see that approaching higher thresholds and redshifts, the slope increases drastically; in particular taking ρ th = 8 when going from z = 0.3 to z = 0.4 the slope is evidently steeper, which results in the much steeper values of the bias in figure 2.
We then plot in figure 4 the magnification bias for third-generation detectors ET and CE.The increased sensitivity of these future experiments ensures that almost all sources are observed, thus the bias will be extremely small.The two detectors differ only towards higher redshifts, with CE showing a steeper tail similar to the LVK case.This could be explained by the different sensitivity curves.CE is sensitive from 5Hz, whereas ET from 1Hz.However, GWs emitted at larger distance will result in a lower merging frequency shown in eq.(2.10), thus a smaller frequency ranges observable by the detector.Therefore, at higher distances,    this might affect the observed number of sources and give rise to the difference seen in the two biases in 4.
For further convenience, we provide functional forms for the magnification biases calculated so far.We fit a polynomial of order 3 with coefficients y = a + bx + cx 2 + dx 3 , and note the fitting values in table 1.Additionally, we opted to fit only the initial part of the curve for LVK, thus before the sudden steepening.As the curves all steepen when the number of observed sources is n obs ∼ 120, we chose this value as upper cut for the functional forms.The appropriate redshift intervals are listed in the final column in table 1.The evolution bias follows similarly.Using its definition in eq. ( 1.2), we can find: We plot b e for a LVK-like experiment in figure 5, and for third-generation detectors ET and CE in figure 6.For the former, the lower sensitivity results in a lower number of detections and thus a much stronger correction.In fact, a strongly negative evolution bias implies that the survey observes a population which appears to become more sparsely distributed as the universe evolves.However, this is actually an effect of the low sensitivity of the experiment, which results in a lower number of detections, and thus it appears as if the sources are intrinsically decreasing with redshift, whilst we should expect a larger volume of sources as we approach cosmic dawn.In fact, the improved sensitivity of 3G detectors ET and CE results in correctly tracing the evolution of BBH mergers.Considering a Madau-Dickinson rate peaking around cosmic dawn as in eq.(2.16), i.e. z c ∼ 2, the evolution bias shows a population of tracers becoming more densely distributed as the universe expands, until it reaches z = 2, where the trend inverses; then, as the tracers start to become more sparsely populated as the universe expands (i.e. the scale factor grows), the evolution bias becomes negative towards present time.Finally, we fit cubic polynomials to the models for b GW s e calculated so far.As with the magnification bias for LVK, we only fit curves in redshifts for which n obs > 120, which we report in the final column of table 2. This was to avoid the final part of the curves, which are drastically impacted by lower observations and thus more uncertain.should be defined in the same manner as for galaxy samples, i.e.: Therefore one only needs to model the number density of SNIa as a function of time and magnitude threshold.

SNIa event rate
In order to model the number density of SNIa, we employ a parameterisation described in [49], constructing a luminosity function Φ(M, z), with M being the peak absolute magnitude of the event.A key assumption in this analysis is that M is given by a narrow Gaussian distribution [50,51].
The starting point of the expression is the star formation rate (SF R).Initially we pick a standard cosmic star formation rate produced [52,53]: Further, the delay time between formation of a binary and subsequent SNIa explosion is modelled by [54] as a power law distribution: Combining the above, we have the SNIa rate in units of yr −1 Mpc −3 [49]: where the factor C SN Ia = 0.032M −1 ⊙ can be computed from the stellar mass range of 3M ⊙ < M < 8M ⊙ for SNIa and the initial mass function [55].Additionally, the explosion efficiency η is taken as the canonical value of 0.04 [52].In the left panel of figure 7 one can see the flat rate density of SNIa as a function of redshift.
Finally, assuming the absolute magnitudes of supernovae are Gaussian-distributed, we can produce the magnitude distribution of sources [49]: where G is a Gaussian distribution, M * = −19.06 in the B-band and σ = 0.56.Φ(M, z) is the magnitude-equivalent of a luminosity function for SNIa.In the right panel of figure 7 one can see the distribution for different redshifts.Here we assume that the dispersion of the distribution is the same remains constant.
The limiting absolute magnitude of detection is found using the limiting apparent magnitude of the detector in a specific band, the redshift of the event and the cross-filter correction K. Recall that at peak brightness: with m and M being, respectively, the apparent and absolute magnitude, and where we assume that the single or cross-filter K-correction is simply a constant offset (note that we neglect errors on the measured apparent magnitude).Thus, the number of sources that have an absolute magnitude smaller than the limiting value -and so are bright enough to be observed -is just the integral of the luminosity function over the allowed magnitudes: Table 3. Coefficients of a third-order polynomial fit to the magnification bias for SNIa surveys with different limiting magnitudes.We report the redshift interval in the final column.Note, the bias is 0 below the redshift range reported for each model.

Magnification and evolution biases of SNIa
As the magnification bias is computed at a fixed redshift, the K-correction can be safely considered constant, thus neglected when differentiating.Thus, the magnification bias can be evaluated by using (3.1): Finally, we find: Thus, using eq.(3.6), we can compute the limiting absolute magnitude at each given redshift from a value of the limiting apparent magnitude of a telescope and compute the magnification bias as prescribed in eq.(3.8) (see figure 8).The evolution bias for SNIa is computed similarly, from eq. (1.2): Both biases are computed up to a value of redshift for which we can observe less than 2σ of the magnitude distribution; this was an arbitrary cut to select only the portion of observations which would be statistically significant.
From figure 8 we note that the magnification bias is fixed to zero until a value of redshift dependent on the magnitude cut.This is due to the fact that before such distance, all SNIa can be observed; however, at such redshift the limiting absolute magnitude approaches the Gaussian distribution of SNIa, thus increasing the number of objects near the threshold (and consequently, of objects being magnified in or out).
As we did for the GWs biases, we fit a cubic polynomial of equation y = a+bx+cx 2 +dx 3 to the values of s SN Ia calculated here, and we report them in table 3.
If we then relax the assumption that SNIa have a fixed brightness across all redshifts [56][57][58][59][60][61] then we can modify (3.6) and add an extra nuisance term ∆m evo (z) on the righthand side, accounting for a potential evolution of the SNIa intrinsic luminosity with redshift. .Left: Magnification biases for a SNIa survey.Solid lines represent the biases for a SNIa with fixed luminosity, whilst dash-dot lines illustrates the bias for SNIa with an intrinsic luminosity redshift dependence.Right Evolution biases for a SNIa survey.Coloured lines show models for an intrinsic SNIa luminosity evolving with time, whilst the black line describe SNIa of fixed peak magnitude M peak = −19.06.In particular, dashdot lines represent an evolution with a power law of index δ 2 = 2, whilst dashed ones have δ 1 = 0.29.For the magnification bias, the two overlap almost completely, whilst this is not the case for the evolution bias, due to extra redshift dependence.Note that the legend is shared between the two plots.
Different models have been proposed to model this evolution, however we will consider only Model B from [57], also illustrated in [62]: We can explore two sets of parameters for this, still in agreement with standard ΛCDM [59], i.e. ϵ 1 = 0.013 ± 0.06 and δ 1 = 0.29 ± 0.22, and ϵ 2 = 0.029 ± 0.052 and δ 2 = 2 ± 1.7.Both imply a population of SNIa which becomes intrinsically dimmer as the redshift increases.The biases are then calculated as before, and then are plotted in figure 8.
If the intrinsic luminosity of SNIa is allowed to vary with redshift, the effect is particularly relevant for the evolution bias, whilst the magnification one is only slightly affected.The reason is due to the nature of the luminosity function dΦ dM in (3.5).For a fixed luminosity kind of SNIa, the magnitude distribution function is separable in redshift and absolute magnitude; hence computing its evolution bias will result in simplifying out the integral over a Gaussian in magnitude (in eq.(3.10)) and only terms related to the intrinsic redshift distribution of SNIa, and its derivative, will remain.These have a small contribution as the redshift distribution is assumed to be relatively flat, as seen in figure 7. Notably, the magnitude cut disappears and thus the evolution bias becomes independent of it.
However, relaxing the assumption of a fixed luminosity brings in an extra redshift derivative of the Gaussian, as now the peak magnitude is dependent on redshift.This allows for the evolution bias to then depend on the magnitude cut as the integral is not simplified out anymore.
Finally, we fit a cubic polynomial to these models for the evolution bias, reporting the coefficients of the functional forms in

Impact on observables
In this section, we will show the relevance of these biases in the luminosity distance clustering power spectra in two ways.Initially, we will investigate their impact on the relativistic corrections to the number counts which are most likely to be detected, namely the Doppler and lensing terms.Then, we will examine how these biases affect the angular power spectrum, exploring both auto-correlation and cross-bin correlations at different redshifts.

Impact on the number counts
The modelling of these biases is extremely important to investigate the relativistic corrections in the number density fluctuation in eq.(1.1).As shown in [37], when analysing clustering in luminosity distance space several terms are dependent on these two parameters, such as the lensing and Doppler terms.Notably, the lensing term in luminosity distance is dependent on both bias parameters, as opposed to the redshift-space case which depends only on the magnification bias.
In particular, by computing the appropriate perturbation in luminosity distance and isolating each correction to the underlying matter density δ n , the number density fluctuation in eq.(1.1) is recast into [37]: where we only report the main correction terms, with coefficients Here r is the source's position, r the comoving distance, γ ≡ rH/(1 + rH), and where we defined Note that β contains the magnification bias s and evolution bias b e of the source type in question.Thus, we explore the impact of the biases we calculated in Sections 2.2 and 3 on the Doppler and lensing terms.We plot the Doppler amplitude A D in figure 9, comparing the results for GWs (for an ET-like experiment) on the left panel and SNIa on the right one.We opt to show only one model of the biases for each population (i.e.Power Law + Peak chirp mass distribution for GWs, and SNIa with fixed intrinsic luminosity) for simplicity.Furthermore, the difference in the Doppler correction with the other models was found to be negligible.Additionally, we plot in grey the Doppler term with, respectively, the biases set to zero (solid line), only magnification (dashed) and only evolution (dash-dot).This clearly shows that the biases dominate the amplitude of the Doppler correction and are crucial for its analysis.In particular, in the case of GWs from an ET-like experiment, the correction traces the evolution bias significantly, whilst for SNIa, the higher values of magnification bias impact the Doppler term and become the main contribution.
We then focus on the lensing amplitude, addressing the same set of biases as in figure 9. On the left of figure 10 we fix the source at z = 2 and explore A L as a function of distance r(z) to the object, thus looking at the integrand on the right-hand side of eq.(4.1).As we previously suggested in [37], the shape of this curve is strictly dependent on the value of the biases, together with the zero-crossing separating magnification near the source from de-magnification away from it.This is clear from the bottom left panel of figure 10 showing the lensing amplitude for SNIa surveys with different limiting apparent magnitudes: higher values of m lim (lower threshold), thus lower magnification bias, shift the zero crossing further from the observer, shrinking the redshift range in which magnification occurs (i.e.A L > 0).We plot the percentage difference between angular power spectra with biases as opposed to without them.The bottom left panel shows models with limiting magnitude m lim > 25 overlapping; bottom right is missing m lim = 25 as this is outside our definition of validity limit for the bias for SNIa (i.e. less than 2σ of the magnitude distribution is observed).Dashed lines represent negative values.

Impact on the angular power spectrum
After examining the impact of the biases on the amplitudes of certain corrections to the number count, we explore their relevance in the angular power spectrum.Using a modified version of the code CAMB4 we produced in earlier work [37], we can compute angular power spectra supplying different values of both magnification and evolution bias.Whilst for GWs we use broad Gaussian windows (σ = 0.2) for the sample bins, motivated by the large uncertainty in the estimation of the luminosity distance of GWs [42,43], for SNIa we chose smaller bins (σ = 0.1), as the related distance uncertainties for LSST will be smaller [63][64][65].
Hence, we plot the percentage difference in the angular power spectrum at each scale when accounting for the biases compared to the case where they are both set to zero for both GWs and SNIa, respectively top and bottom of figure 11.We do this for three different redshifts, noting that z = 1.5 is outside what we defined as the validity range of the bias for SNIa assuming a limiting apparent magnitude m lim = 25, and thus this particular case is not shown in the bottom right panel of figure 11.We recall that the SNIa validity limit was fixed at the distance at which 2σ of the magnitude distribution is observed.
One can immediately note the difference between the two tracers.Whilst SNIa show an increasing strength of the biases at higher redshifts, the models considered here for GWs impact differently.In fact, at low redshift the effect of the biases is simply at a percentage level, approaching a negligible one around z = 1.5, and rising steeply after.This follows directly from the models of the biases studied here: the evolution bias (see figure 6) is close to −2 at very low redshifts, crosses 0 around z = 1.5, and then rises further.
We then examine the impact of the models of the biases considered when cross-correlating different redshift bins.We show two different examples of this: setting a background tracer at z b = 1.5 and a foreground one at z f = 0.5, and similarly with z b = 2.0 and z f = 1.0.This is shown in figure 12, with GWs in the top panels and SNIa in the bottom ones.Similarly to figure 11, the plots report the percentage difference between the angular power spectra with biases compared to those with the biases set to zero.It is clear how not accounting for the magnification and evolution biases can lead to large differences, depending on the tracer considered.On the top left panel, cross-bins angular power spectra for GWs show percentage level difference, although the contrast increases drastically by roughly an order of magnitude when shifting to higher redshifts (on the top right).On the other hand, angular power spectra built using SNIa already show differences of several times those with biases set to zero; going to higher redshifts can push the percentage difference to a few order of magnitude.In all cases, the difference tends to a constant value at a larger value of ℓ which depends on the redshift considered.
The large differences clearly show that cross-correlations require considering the impact of the magnification and evolution biases for both types of transient tracers.Additionally, cross-bin correlations allow for the different bias models to be distinguishable even at lower redshifts.This is clear from the plots on the left-hand side in figure 12, where each model produces a distinct difference with respect to the angular power spectrum without biases.Such contrast is not seen in the auto-correlations with GWs in the top panels of figure 11 even at z = 1.5, where the two models are both roughly of the same order of magnitude.SNIa models already showed substantial differences in auto-correlations at high redshift, as shown in the bottom right panel of figure 11.However, the distinction becomes more significant in the cross-bins correlation, with differences of even an order of magnitude arise between different limiting magnitude models, as shown in the bottom plots of figure 12.
Finally, we explore the impact of the biases modelled in this paper on the auto-correlation angular power spectrum across redshift.We plot the percentage difference between the nonzero biases case with respect to the zero bias case as a function of redshift in figure 13.The step in the plot represents the binning used for both GWs and SNIa, respectively z = 0.2 and z = 1; thus, we decided to avoid a smooth curve for clarity.On the LHS we plot the difference when calculated for GWs observed with an ET-like detector, and on the RHS the same for SNIa from a magnitude-limited survey.The former shows (positive) deviation from zero from around z = 1, growing almost linearly with redshift.It also shows that before such a point, the angular power spectrum with biases accounted for grows slightly and turns back towards zero.An explanation of this could be traced to the behaviour of the magnification and evolution biases for ET and CE (shown in figure 6).Since s GW is significantly close to zero, the largest contribution will be given by the evolution bias; and considering that b GW e crosses zero around z < 1.5, this might explain the dip in the percentage difference for GWs in figure 13.
On the other hand, as shown in figure 8, the evolution bias of SNIa was very small, whereas s SN Ia was substantially larger than the GWs case.The Doppler kernel in figure 9 and the lensing amplitude in figure 10, the correction terms for SNIa were dominated by the impact of the large magnification bias; similarly, the difference in the auto-correlation angular power spectra between the case with biases and the one with them set to zero, is significantly   As before, we plot the percentage difference between the C ℓ with biases as opposed to those without.
impacted by the magnification bias.This is clear by the shape of the curves shown in RHS of figure 13: the difference is zero until the redshift value at which the corresponding s SN Ia starts growing.
Figure 13 further stresses that for SNIa the models can be clearly distinguished from z ∼ 1; however, whereas auto-correlation angular power spectra with GWs do highlight the need for accounting for the biases, they do not provide large enough differences between the models considered.Such distinction is only seen when cross-correlating separate redshift bins, as seen in figure 12.

Summary and conclusions
By modelling and ultimately measuring the bias properties of tracers, we can use them as probes of large-scale structure and cosmology.The bias properties of galaxies is a well-studied topic; with the advent of large data sets of transient sources, we need to learn how to describe their biases equally well.As we have eluded to, there are particular fundamental physics motivations for wanting to use tracers who most naturally 'belong' in luminosity distance space, such as GW sources and SN1a.
The observed fluctuation of the source number counts, ∆ O , is the basis of many clustering analyses [7,8].When considering transient tracers such as GWs or SNIa, this has to be computed in luminosity distance space, as the two types of object do not carry direct information of their redshift [37].The expression for the number counts contains several relativistic corrections, notably Doppler magnification and lensing.These depend on the magnification and evolution biases, which can significantly alter their amplitudes.
In this paper, we first described the modelling of magnification and evolution biases of gravitational waves from binary black hole mergers.We initially carefully described the necessity of defining them with respect to a specific detector and not to any single "intrinsic" quantity, as opposed to the traditional case in galaxy clustering.We also employed two different distributions of chirp masses, namely a Power Law + Peak and a Broken Power Law consistent with population analysis of the GWTC-3 catalogue by the LVK collaboration [42,43].Furthermore, we examined GW biases for three different SNR thresholds, noting that a higher ρ th results in larger values of the biases.The lower sensitivities of the present terrestrial detectors (relative to 3G expected sensitivity curves) yields a strongly negative evolution bias.This implies that sources are becoming sparser with redshift; however, given the expectation of a peak around cosmic dawn (z ∼ 2), it is clear this is due to a decline in the detector's sensitivity at these redshifts.In fact, third-generation observatories ET and CE should instead be able to perfectly trace the evolution of GWs from BBH mergers.Consequently, their related magnification biases are particularly small up to high redshift.We note that the method described for GWs from BBH mergers can be equally applied to different sources of GWs, such as neutron stars mergers or black hole-neutron star binaries, provided the appropriate chirp mass distribution and merger rate are supplied.Further, we explored the biases for SNIa, investigating the impact of changing the magnitude threshold.Whilst the magnification bias is especially sensitive to this, the evolution bias for the SNIa models considered is independent of it.Additionally, the former can reach much larger values than the ones seen for GWs, while the latter is close to flat around zero, implying a population density that remains close to constant across redshift.We also explored the possibility of a population of supernovae of intrinsic luminosity evolving with redshift, using two different models proposed in the literature.We found that the magnification bias is only slightly altered, while the evolution one changes significantly depending on the model adopted (i.e.whether the intrinsic luminosity of SNIa decreases or increases with redshift).
Finally, we investigated the effects of these biases on relativistic corrections to the observed number counts and on the angular power spectrum.In both cases, the impact of the biases produces a significant difference with respect to cases where the biases were set to zero.We found that the lensing magnification is strongly dependent on the magnification bias, as expected.However, the relativistic Doppler correction shows a different behaviour depending on the tracer.For GWs it is sensitive to the evolution bias, whilst for SNIa the dominant bias is the magnification bias.This is explained by the different values taken by both s and b e for the two different tracers: the former being significantly small for GWs, the latter for SNIa.Therefore, with one of the two parameters closer to zero, the other has a stronger impact.
The same effect occurs then when comparing angular power spectra with biases included to ones with biases set to zero, as in figure 13.This in turn would have an impact on any derived constraints from the angular power spectra in luminosity distance.For GWs, the difference between the angular power spectra is initially small, rising after the redshift for which b GW s e > 0. For SNIa, the difference is null until the distance at which s SN Ia > 0. Furthermore, when analysing the impact of the biases on the angular power spectrum, we found a percentage level difference in auto-correlations with GWs between spectra with biases accounted for and spectra with biases set to zero.When the same analysis is applied to SNIa, the differences are greatly increased at high redshifts, i.e. when s SN Ia > 0. This shows the importance of accounting for the biases when computing angular power spectra, however, at least for GWs, the different models of the biases are distinguishable only at high redshifts.The separation between models is instead achieved more clearly when investigating cross-bin correlations between different redshifts.Figure 12 not only shows greatly increased differences with respect to the unbiased spectra, but also highlights the possibility of distinguishing the specific model of the bias.As these parameters are strictly dependent on population properties, distinguishing between them could help us infer details of the mass distribution of these tracers.
Having set up frameworks with which to model transient biases, our next step will be to understand how these impact constraints from cross-correlations in the era of stage IV galaxy surveys, and the 3G era of gravitational wave detection.This impacts not only 3G cosmology, but may also be highly relevant for astrophysics and compact object formation channels, if we find that population properties can be simultaneously constrained.no.URF\R1\180009).CC is supported by the UK Science & Technology Facilities Council Consolidated Grant ST/T000341/1.

B Chirp Mass distribution
We want to compute the PDF h(z) of the chirp mass M given the distributions of primary and secondary masses, g(m 1 ) and f (m 2 ).Let us recast these random variables (for simplicity) in m 2 = x,m 1 = y, M = z, and their relation is given by the chirp mass definition: Now, looking at (B.1), we need to find the solutions for x(z, u).We can recast the equation as a third order equation for x of the form

Figure 1 .
Figure1.Left: Fraction of sources detected, S(ρ th /ρ 0 ) as a function of ρ 0 , i.e. the characteristic SNR of the source.Low values of S imply either a low SNR threshold or loud events, whereas large S signifies a high threshold or quiet events.Center: Fraction of sources undetected at each random orientation.θ is equivalent to the ratio between the SNR threshold ρ th and the characteristic SNR of the source, ρ 0 .Right: Orientation PDF.

Figure 2 .
Figure2.Magnification bias for an LVK-era experiment.We stop calculating s for each model when the number of observed sources is below the arbitrary value of 100 events Gpc −3 ; this is to suggest a limit of validity for the bias, as clustering analysis requires a larger number of sources.Solid lines indicate a Power Law + Peak model for the primary BH mass, whilst dashed ones show the Broken Power Law model.Fainter lines represent the range where the biases become more uncertain, as the number of observed sources crosses 120, getting closer to our arbitrary cut.

Figure 3 .
Figure 3. Modelling the number densities of GWs observable with LVK using (2.5).Left: as a function of redshift for different (fixed) SNR threshold values; Right: as a function of ρ th at different (fixed) values of redshift.At each given z, the value of the magnification bias is given by the slope of the corresponding line in the right-hand plot.Increasing ρ th starts to change significantly around z = 0.5, explaining the behaviour of s in figure 2.

Figure 4 .
Figure 4. Magnification bias for third-generation detectors, left ET, right CE.As before, solid lines are for a Power Law + Peak distribution of the primary BH mass, and dashed ones are for the Broken Power Law model.

Figure 5 .
Figure 5. Evolution bias for a LVK-like survey.As previously, solid lines indicate a Power Law + Peak model for the primary BH mass, dashed for the Broken Power Law model.Similarly to figure 2, faint lines represent a regime with observed number of sources 100 < n obs < 120, thus approaching our arbitrary cut for the biases.

Table 2 .
04 −1.76 × 10 1 1.05 × 10 2 −4.36 × 10 2 [0.1, 0.22] Coefficients of a third-order polynomial fit to the evolution bias for GWs detectors for different SNR threshold.Given the negligible difference in the biases between ET and CE over the redshift interval z ∈ [0.1, 3.5], we only report the fits to ET.

Figure 6 .Figure 7 .
Figure 6.Evolution bias for third-generation ground based GWs observatories ET (left) and CE (right).These will be able to trace the evolution of BBH mergers up to high redshift.

27 Figure 8
Figure 8. Left: Magnification biases for a SNIa survey.Solid lines represent the biases for a SNIa with fixed luminosity, whilst dash-dot lines illustrates the bias for SNIa with an intrinsic luminosity redshift dependence.Right Evolution biases for a SNIa survey.Coloured lines show models for an intrinsic SNIa luminosity evolving with time, whilst the black line describe SNIa of fixed peak magnitude M peak = −19.06.In particular, dashdot lines represent an evolution with a power law of index δ 2 = 2, whilst dashed ones have δ 1 = 0.29.For the magnification bias, the two overlap almost completely, whilst this is not the case for the evolution bias, due to extra redshift dependence.Note that the legend is shared between the two plots.

. 4 )= 27 Figure 9 .
Figure 9. Dimensionless Doppler kernel A D including the contribution of the biases.Left: the effect when considering GWs observed by ET (given the similarities of the biases between ET and CE we plot only the former); Right: the effects when considering SNIa with different magnitude cuts.For both cases, we plot in grey the contribution in the case of the biases set to zero (solid line), only non-zero magnification bias (dashed) and only non-zero evolution bias (dash-dot); in particular, the latter two are set for ρ th = 8 on the left, and m lim = 26 on the right.

Figure 11 .
Figure 11.Impact of the models of the biases considered on auto-correlation angular power spectra for GWs as seen by an ET-like experiment (top) and for a magnitude limited SNIa survey (bottom).We plot the percentage difference between angular power spectra with biases as opposed to without them.The bottom left panel shows models with limiting magnitude m lim > 25 overlapping; bottom right is missing m lim = 25 as this is outside our definition of validity limit for the bias for SNIa (i.e. less than 2σ of the magnitude distribution is observed).Dashed lines represent negative values.

27 Figure 12 .
Figure 12.Impact of the models of the biases considered on cross-bins correlations, plotting the percentage difference between angular power spectra with biases and those without.Left: correlating a background bin at z b = 1.5 with a foreground one at z f = 0.5; right: the same with background one set to z b = 2 and foreground one to z f = 1.As in figure11, top is for GWs for an ET-like experiment, and bottom for SNIa.The bottom plots lack the model with limiting magnitude m lim = 25 as it is outside the validity limit.

Figure 13 .
Figure 13.Impact of the biases considered on the angular power spectrum at ℓ = 10 across redshift.As before, we plot the percentage difference between the C ℓ with biases as opposed to those without.

2 )
The joint PDF p(z, u) of the new variables z, u is related to that of x, y by the conservation of the volume in the space of probability f (x)g(y)dxdy = p(z, u)dudz ,(B.3)where we have assumed that f (x), g(y) are independent.The differentials are related by the determinant of the Jacobian matrix of the transformation:dxdy = |J(z, u)|dudz .(B.4)Taking x i as a root of the equation z = z(x, y), and the fact that, from (B.2), y(z, u) = u, we can write the Jacobian as:|J(z, u)| = ∂x i (z, u) ∂z .(B.5)Substituting (B.5) and (B.4) into (B.3), and summing over the possible roots, we find the combined PDFp(z, u) = x i f (x i (z, u))g(u) ∂x i (z, u) ∂z , (B.6)and to obtain the PDF for z, h(z), we marginalise over u: h(z) = x i du f (x i (z, u))g(u) ∂x i (z, u) ∂z .(B.7)

x 3 + 5 u 3 , c = − z 5 u 2 . (B. 8 )
bx + c = 0 , with b = − z Depending on the sign of the delta of the cubic equation above

Table 1 .
Coefficients of a third-order polynomial fit to the magnification bias for GWs detectors for different SNR threshold.

Table 4 .
Coefficients of a third-order polynomial fit to the evolution bias for SNIa surveys with different limiting magnitudes.As before, We report the mean squared error (MSE) and redshift interval in the final two columns.Note, the first model refers to a population of SNIa of fixed intrinsic luminosity, and yields a result independent of the limiting magnitude.Models δ 1 and δ 2 describe SNIa with intrinsic luminosity evolving with redshift at different rates.