{'fraction_non_alphanumeric': 0.02768795472918351, 'fraction_numerical': 0.01864389652384802, 'mean_word_length': 4.476618705035971, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The redshift dependence of black hole mass distribution: Is it reliable for standard sirens cosmology?

An upper limit on the mass of a black hole set by the pair-instability supernovae (PISN) process can be useful in inferring the redshift of the gravitational wave (GW) sources by lifting the degeneracy between mass and redshift. However, for this technique to work, it is essential that the PISN mass-scale is redshift independent or at least has a predictable redshift dependence. We show that the observed PISN mass-scale can get smeared and the position of the PISN mass-scale is likely to exhibit a strong redshift dependence due to a combined effect from the non-zero value of the delay time between the formation of a star and the merging of two black holes and the metallicity dependence of PISN mass scale. Due to the unknown form of the delay-time distribution, the redshift dependence of the PISN mass cut-off of the binary black holes (BBHs) cannot be well characterized and will exhibit a large variation with the change in redshift. As a result, the use of a fixed PISN mass scale to infer the redshift of the BBHs from the observed masses will be systematically biased. Though this uncertainty is not severe for the third observation run conducted by the LIGO-Virgo-KAGRA collaboration, in the future this uncertainty will cause a systematic error in the redshift inferred from the PISN mass scale. The corresponding systematic error will be a bottleneck in achieving a few percent precision measurements of the cosmological parameters using this method in the future.


INTRODUCTION
Inference of the cosmological parameters from the gravitational wave (GW) sources is one of the key science goals of the currently ongoing network of GW detectors (Abbott et al. 2018) such as LIGO (Aasi et al. 2015), Virgo (Acernese et al. 2014), and for the upcoming GW detectors such as KAGRA (Akutsu et al. 2020), LIGO-India (Unnikrishnan 2013), LISA (Amaro-Seoane et al. 2017), Cosmic Explorer (Reitze et al. 2019;Hall & Evans 2019), and Einstein Telescope (Punturo et al. 2010), as it can provide accurate measurement to the luminosity distance of the GW sources within the framework of the general theory of relativity, and without invoking any additional distance calibration. This was shown for the first time in the seminal work by Schutz (1986), which justifies the reason for calling GW sources the standard sirens. However one of the essential requirements for making robust measurements of the cosmological parameters from standard sirens is to be able to make an independent and accurate inference of the redshift of these sources. Standard sirens, though excellent distance tracers, cannot provide the redshifts to the sources independently, unless one can break the degeneracy between mass and redshift, by using a known mass scale. For binary neutron star (BNS) sources, one can use the tidal † smukherjee1@perimeterinstitute.ca deformability to break the mass-redshift degeneracy (Taylor et al. 2012;Messenger & Read 2012). On the other hand for the binary black holes (BBHs), the existence of a maximum mass of the black holes (BHs), can be used to infer the redshift to the sources You et al. 2021;Mastrogiovanni et al. 2021a). From the theory of stellar models of BH formations, it is predicted that there is a maximum mass of the BHs due to the pair-instability supernovae process (PISN) (Bethe & Brown 1998;Portegies Zwart & Yungelson 1998;Belczynski et al. 2002), and the maximum mass of the BH is expected to be between 40 − 50 M Renzo et al. 2020). A recent study has also shown , that the variation of the maximum mass-scale changes by only about 10% due to a large variation in the stellar metallicity and stellar winds. As a result, if the PISN mass-scale is robust within 10% accuracy, it can be used to break the mass-redshift degeneracy and can be used to infer the cosmological parameters which affect the cosmic expansion history. Even in the absence of an accurate theoretical prediction of the maximum mass value, one can do a joint estimation of the cosmological parameters and the PISN mass-scale You et al. 2021;Mastrogiovanni et al. 2021a). So, this particular method can be extremely powerful also for the third-generation GW detectors such as Cosmic Explorer (Reitze et al. 2019;Hall & Evans 2019), and Einstein Telescope (Punturo et al. 2010) to infer the cosmic expansion history to a high redshift z ≈ 50, which is unlikely from any electromagnetic probes in the same time scale.
However, for the usage of the PISN mass-scale for the cosmological purpose, the mass-scale must be redshift independent, or at least the redshift dependence should be predictable. If the PISN mass scale evolves with cosmological redshift, then our ignorance of the redshift dependence on the PISN mass scale, will lead to a systematic error in the measurement of the true cosmological redshift and hence will affect the inference of cosmological redshift that will use this redshift. Several previous studies assumed that the PISN mass-scale is redshift independent or will exhibit mild redshift dependence based on the theoretical studies Renzo et al. 2020), and as the stellar metallicity at low redshift (z < 2) varies not very significantly. In this work, we scrutinize whether the assumption that the PISN mass-scale is redshift independent is valid and show whether we can use it reliably to infer the true redshift to the binary BHs. Recent studies have explored the variation of the mass distribution with redshift for a fixed model of cosmology from GWTC-2 (Fishbach et al. 2021) and GWTC-3 (Abbott et al. 2021a) and did not find any statistically significant deviation from the redshift independent scenario.
We show that even though the variation of the PISN massscale is below 10% across a large variation in the metallicity range (Z ∈ {10 −3 − 10 −5 }) , the delay time between the formation of the stars and merger of the BBHs will play a key role in the observed mass distribution of the BBHs. The observed mass distribution of the BBHs at a redshift will arise from the BHs formed over a wide range of cosmological redshifts due to a non-zero value of the delay time between the formation of star and merger of BHs. As a result, the mass distribution of the BBHs at any redshift will come from a vast range of redshifts over which the stellar metallicity can vary significantly. In this paper, we show the impact of the observed mass distribution of BBHs on the inferred redshift if a fixed value of the PISN mass scale MPISN is assumed.
The paper is organized as follows, in Sec. 2 we introduce the idea behind the redshift dependence of the PISN mass scale and show its impact on the GW source population in Sec. 3. In Sec. 4 we show the inference of the PISN mass scale from the mock GW samples for the network of LIGO-Virgo-KAGRA (LVK) detectors and Cosmic Explorer (CE) and in Sec. 5 we discuss its impact on the GW population and the cosmological parameters. Finally in Sec. 6, we discuss the conclusion and the future prospects.

LOWER EDGE OF THE PISN MASS-GAP AND ITS REDSHIFT DEPENDENCE
The gap in the mass distribution of stellar-mass BHs is expected due to the mass loss of the heavy stars due to the PISN process (Fowler & Hoyle 1964;Rakavy & Shaviv 1967;Fraley 1968;Bond et al. 1984;Woosley et al. 2002;Talbot & Thrane 2018;Marchant et al. 2019;Stevenson et al. 2019). The lower limit of the mass gap is expected to be around 45 M , which is set by the mass loss during the PISN process Renzo et al. 2020). The value of the lower limit of the PISN mass is shown to vary by only about 7% for the variation in the stellar metallicity Z from 10 −5 to 3 × 10 −3 , as the mechanism of wind-loss depends strongly on the stellar metallicity (Vink et al. 2001;Mokiem et al. 2007).
Though multiple other parameters such as the nuclear rates can affect the PISN mass-scale MPISN = 45 M by about 35%, this parameter is not expected to evolve with cosmological redshift and hence cannot cause redshift evolution of the PISN mass-scale. So, the stability of the lower edge of the PISN mass gap can be a robust feature for cosmology. However, whether the PISN mass-scale can be used for reliably inferring the cosmological parameter depends on whether the observed mass distribution is redshift independent or not, or at least whether the redshift dependence of the PISN mass-scale can be well characterized. In this work, we scrutinize the possible redshift dependence of the observed mass distribution of BBHs even in the scenario when the PISN mass scale is well predicted by theoretical studies. Though in reality, there may be uncertainties in the theoretical understanding as well. As a result, the uncertainty that we are exploring in this paper should be considered as only a lower bound. There can be additional errors on top of this as well.
We model the metallicity dependence of the PISN mass cut-off by a relation where a weak dependence on the metallicity can be captured as a logarithmic correction with a free parameter α. For α = 1.5 at Z * = 10 −4 and MPISN(Z * ) = 45 M , we can capture the variation in the PISN mass scale shown in  due to the change in metallicity over the range Z ∈ [10 −5 , 5 × 10 −3 ]. Along with the evolution of the metallicity, the PISN mass scale can also vary due to changes in the fundamental physics such as the reaction rates (Mehta et al. 2022). The stellar metallicity in the Universe evolves with redshift (Mannucci et al. 2010;Sommariva et al. 2012;Krumholz & Dekel 2012;Dayal et al. 2013;Madau & Dickinson 2014). The metallicity at a high redshift (z > 2) is much smaller in comparison to the low redshift Universe z < 2. The firstgeneration stars contaminate the interstellar medium and cause a chemical evolution of the Universe. We can treat the metallicity evolution with redshift by a relation log 10 (Z(z)) = γz + ζ, where γ captures the redshift dependence and ζ captures the metallicity value at z = 0 (Mannucci et al. 2010;Madau & Dickinson 2014). This relation captures the metallicity of the parent star or the gas cloud from which a star has formed. It is written to express only a mean evolution of the metallicity. Along with the mean metallicity evolution of the Universe, there is going to be a scatter in the metallicity depending on the galaxy properties. Such a source of uncertainty brings additional stochasticity to the metallicity relation. Currently, a limited number of observations (Gallazzi et al. 2008;Mannucci et al. 2010;Krumholz & Dekel 2012) are available to explore the environment dependence of the metallicity, and most of our current understandings are based on simulations (Genel 2016;Torrey et al. 2019). These studies show that the overall median metallicity dependence of the galaxies at different can be explained by power form (Pei et al. 1999  studies of GW merger rates and mass distribution are performed (Belczynski et al. 2002;Dominik et al. 2012;Dominik et al. 2015;Mapelli et al. 2017;Giacobbo et al. 2018;Toffano et al. 2019;van Son et al. 2022) which are motivated by these studies and show the black hole mass distribution can exhibit a redshift dependence. The existence of any stochasticity in the galaxy metallicity distribution will also influence the mass distribution but is currently not well known. However, as the relation given in Eq.
(1) is in terms of the logarithm of metallicity, so the impact of fluctuation around the median value depending on the individual galaxy properties is going a small (logarithmic) change. As we are unable to measure the host of the BBH due to a large sky localization error of the BBH, we cannot directly associate the properties of galaxies with BBH source properties. So, we can only infer an ensemble average mass distribution from the GW data and the additional stochasticity (which will depend on the host properties) will appear as an additional uncertainty in the measurement of MPISN. As a result, we consider a median distribution of galaxy metallicity and the dependence of MPISN on it.
The observed mass distribution of the BBHs at a redshift zm is going to contribute from the BHs which have formed at an earlier redshift z < zm due to a non-zero value of the delay time t d . The formation time of the individual companion objects will be at a different redshift. As a result, the metallicity dependence of the star corresponding to that will play a role. So, for a probability distribution P (t d ) of the delay time, there is going to be mixing between the BHs forming at different redshifts. The corresponding window function of masses of the BBHs mergers can be written as where N is a factor to normalize the window function. The term Pt d (zm, z ) denotes the probability for a BBH source to the merger at a redshift zm formed at a redshift z with a delay time t d , and the term Ws(m(z )) denotes the window function of the BH masses at a redshift z with the mass m denoting the mass of the individual objects in their source frame. The convolution of the probability distribution of the delay time and the probability distribution of the mass distribution gives the source frame mass distribution at redshift zm.
The mass distribution of BBHs at redshift z can be modified depending on the metallicity evolution in the Universe. We show a schematic diagram in Fig. 1 explaining the variation in the PISN mass of the BHs merging at a fixed redshift. At a fixed redshift (denoted by z), the merging BBHs could have formed at different redshifts (z1', z2', z3', and z4') from the parent stars at redshifts (denoted by redshifts z1, z2, z3, and z4) 1 . So, the PISN mass scale of the BHs formed at different 1 The time lapsed between the formation of a star and the BH can be negligible in comparison to the cosmic time scale of it to merge (which is a few hundred Myr to Gyr) for the heavier black holes that contribute to the PISN mass scale. The redshift of the two parent stars are nearly the same z2 ∼ z1 (and also z4 ∼ z3). However, the difference between the redshifts like z3−z1 of individual binary pairs can be large.
redshifts will be different for the heavier BH. The probability distribution of the delay time with a non-zero minimum value leads to the mixing of BHs formed at a different cosmic time to merge with different PISN masses. As a result, BBHs will not have a fixed PISN mass scale and can show a large variation even if the variation due to a change in the metallicity of the PISN mass is small. The signature of the delay time distribution on the BHs mass distribution can arise in all those scenarios of binary formation for which the delay time distribution can go beyond a few hundreds of Myrs. Only for the scenarios where BBHs are merging very fast (less than a few tens on Myrs), the masses of the BBHs may not show a significant difference, if the metallicity of host galaxies are similar. The delay-time distribution and the corresponding value of the minimum delay time are not well known. Several theoretical studies are made to understand the delay time distribution and its relation with the BBH formation channels (O'Shaughnessy et al. 2010;Banerjee et al. 2010;Dominik et al. 2012;Dominik et al. 2015;Mandel & de Mink 2016;Lamberts et al. 2016;Cao et al. 2018;Vitale et al. 2019;Elbert et al. 2018;Eldridge et al. 2019;Callister et al. 2020a;du Buisson et al. 2020;Santoliquido et al. 2021;Safarzadeh et al. 2020;van Son et al. 2022). For a uniform in the log-space distribution of the separation between the binaries, there is going to be a power-law form of the delay time distribution denotes the minimum value of the delay time. For t d < t min d , the probability distribution is zero. However, studies have shown that the probability distribution of the delay time can have a different functional form, apart from the unknown value of the minimum delay time denoted by t min d . Recently, constraints on the delay time are obtained from the individual events (Fishbach & Kalogera 2021) and the stochastic GW background (Mukherjee & Silk 2021). In the future, it will be possible to measure this quantity to much higher precision by correlating the astrophysical properties of the GW sources with the star formation and metallicity properties of the host galaxies, which can be inferred from different emission lines signatures (Mukherjee & Dizgah 2021). Gravitational lensing of GW can also provide an independent way to measure the delay time distribution from the current generation detectors (Mukherjee et al. 2021e).
For large delay time scenarios, the component mass of the BHs can arise from a very high redshift. At high redshifts, the PISN mass-scale is going to be different and those sources contribute to the BBHs merger. For a probability distribution of the delay time P (t d ) = t −κ d for t d t min d , we will witness a mixing of the BHs from different redshifts to contribute at the redshift zm. The probability mass distribution is going to depend on parameters κ, t min d , γ, α. We show the probability distribution of mass at different redshifts for a few different values of the parameters such as κ, t min d , γ, α in Fig. 2. For the values chosen for the parameters κ, t min d , γ, and α, the variation in the PISN mass scale is well within the variation in the mass range explored previously (Belczynski et al. 2002;Belczynski et al. 2016;Mapelli et al. 2017;Giacobbo et al. 2018;Toffano et al. 2019;Farmer et al. 2019). Fig. 2 shows that for different redshifts, the window function for the GW mass-scale is not the same at all redshifts. At higher redshifts, the PISN mass-scale exhibits a larger value of the PISN mass cut-off than at a lower redshift. This vari-ation in the mass scale is more if the metallicity evolution with redshift is large or the delay time distribution is large or with a shallow evolution of the probability distribution (κ close to zero). The shift in the PISN mass scale MPISN and the maximum mass Mmax with redshift is shown in Fig. 3. It indicates that for different scenarios, there is a significant variation in the mass scales. This implies that the observed mass distribution of the BBHs is going to be significantly different because of the variation of the parameters related to metallicity and delay time distribution. For the scenarios considered in this analysis, the variation in PISN mass scale MPISN can be between 10% to 50%. The variation in the maximum mass scale Mmax can be also around 25%. A previous study explored the metallicity dependence of the Mmax (Safarzadeh & Farr 2019) but did not show the impact of delay time on the mass distribution and how it can naturally produce a redshift dependence of the PISN mass scale and the smearing property of the BH mass distribution. The corresponding source frame probability mass distribution merging at a redshift zm can be written as P(m(zm)) = W (m(zm))Ps(m(zm)), where Ps(m(zm)) is the probability distribution of the stellar compact objects that source the BH formation. This probability distribution can be considered to be motivated by a simple power-law form m −α . We assume in this paper that Ps(m(zm)) is redshift independent. But a redshift-dependent scenario is not outside the theoretical scope. The observed masses of the GW sources merging at redshift zm will be redshifted which can be written as We show the probability distribution of the masses, for this simple power-law form in Fig. 4 for two different redshift values z = 1.1 and z = 6.1. So, even though the mass distribution of the BHs is extended to higher masses, the mass window function W(m(zm)) suppresses the masses above the mass scales set by the window functions at those two redshifts. The mass spectrum shows that the cutoff happens at different mass values at different redshifts. As result, the observed mass distribution of BHs is no more redshift independent.
In summary, even a marginal dependence of the PISN mass scale of a BHs can lead to a much broad distribution of the PISN mass scale due to the mixing between the BHs formed at a different cosmic time due to a non-zero value of the delay time between the formation of a star and merger of a black hole. The probability distribution of the delay time leads to redshift dependence in the observed PISN mass scale and the maximum mass of a BH. We summarise below the key aspects of the redshift dependence of the mass distribution.
The key aspects of this mass modeling • This model predicts one of the potential sources of redshift dependence of the PISN mass scale MPISN and shows how it is correlated with the quantities such as stellar metallicity, delay time distribution, and the value of the minimum delay time.
• This model naturally predicts a smearing behavior between the PISN mass scale and the maximum value of the mass distribution and how it evolves with redshift.
• This model predicts a redshift dependence of the maximum mass value denoted by Mmax and also its correlation In light of the recent results from the LVK observations Abbott et al. (2019Abbott et al. ( , 2020Abbott et al. ( , 2021b, we have detected a few events whose masses are higher than the usual PISN scale considered as 45 M . As we have shown in Fig. 2, due to a non-zero value of the delay time and evolution of the stellar metallicity, the mass distribution of the BHs can go up to high redshift. The window function W(m) gets a slope in its distribution, instead of a sharp cutoff. The slope of the mass window function can go beyond the value 45 M , which can allow for heavier masses of the binary mergers. High mass systems can be also be possible from several astrophysical scenarios as shown previously (Di Carlo et al. 2020;Farrell et al. 2021;Costa et al. 2022;Briel et al. 2022). It is also worth mentioning that apart from the first-generation BHs, second-generation BHs due to hierarchical mergers can lead to heavier BHs. This can lead to additional variation in the BH mass distribution (Liu & Lai 2021;Barrera & Bartos 2022). In summary, all these astrophysical uncertainties can make the mass distribution of the BHs redshift dependent which can make it more difficult to use for the cosmological purpose.

IMPACT ON THE MASS DISTRIBUTION OF THE GW SOURCES
The redshift evolution of the PISN mass scale with redshift will lead to a change in the mass distribution of BBHs with redshift. We study the impact on the observed mass distribution of the binary sources for two different network configurations, the current generation detectors such as LIGO (Aasi et al. 2015), Virgo (Acernese et al. 2014), and the upcoming GW detectors such as KAGRA (Akutsu et al. 2020) and LIGO-India (Unnikrishnan 2013), and the third generation GW detectors such as Cosmic Explorer (Reitze et al.  (Punturo et al. 2010). The sources detectable from the current generation detectors and the next generation detectors will have a redshift reach up to z = 1 and z ∼ 80 respectively. In our analysis, we primarily focus on the GW sources for redshift up to z = 7 for the third-generation detectors, since the expected merger rate at higher redshifts is not well known and difficult to model for the astrophysical BHs which are motivated by the Madau-Dickinson star formation rate (Madau & Dickinson 2014). We sample GW sources with a mass distribution of the power-law form m −α with α = 2.35 and a cut-off at a mass scale according to the mass window function which is a function of redshift. Along with a power-law form, we consider a bump at the MPISN scale with the relative height of the peak in comparison with the low mass value as Λg = 0.1, which is in the agreement with the GWTC-3 observation from LVK collaboration (Abbott et al. 2021d,a,b,c) and also with the previous population synthesis models (Stevenson et al. 2019). The sources at higher redshifts have a mass cutoff at a higher value resulting in a broader mass spectrum at a high redshift than at a low redshift. We use a network of four GW detectors (LIGO-Hanford, LIGO-Livingston, Virgo, and KAGRA) to estimate the posteriors on the mass samples for GW sources distributed up to redshift z = 1.1 which is possible to access from the current generation detectors (Abbott et al. 2018).

2019; Hall & Evans 2019), and Einstein Telescope
We use an approximated form of the likelihood to infer the posteriors on the GW source parameter, following the previous analysis Mastrogiovanni et al. 2021a). This assumption will not change the main conclusion of the paper on the incorrect estimation of the redshift using the PISN mass scale. The matched filtering signal to noise ratio (SNR), ρ is estimated for BBHs with detector-frame chirp mass M d at the luminosity distance of d , using the relation (Farr et al. 2016;Mastrogiovanni et al. 2021b) where Θ is the detector projection factor that we assume as uniform between [0, 1]. We have made a pessimistic choice of the projection factor to show the impact of the bias in inferring the redshift using GW sources. For sources, with a higher value of the projection factor, the bias in the redshift estimation can be even more noticeable. The value of ρ * = 8 is set for the parameters M * d and d * for the LVK design sensitivity (Abbott et al. 2018) and CE (Hall & Evans 2019). The network matched filtering SNR is obtained by using ρ 2 det = i ρ 2 i . We choose BBHs having ρ det 10 as detected events. The estimated mass distribution as a function of redshift shows the variation due to the evolving mass distribution. This indicates that the assembling of the heavier objects takes place at higher masses for higher redshifts. The corresponding redshift evolution of the mass distribution is shown in Fig. 5 only for the events which are close to the PISN mass scale.

INFERENCE OF THE SOURCE REDSHIFT USING THE PISN MASS CUTOFF
The PISN mass scale can be used to identify the redshift of the GW sources if the mass of the PISN mass scale is redshift independent. However, as we have discussed in the previous sections, the PISN mass scale is going to vary with redshift. As a result, if we would like to infer the redshift of the GW sources for a fixed value of the PISN mass scale, then there is going to be a systematic error in the redshift estimation. The redshift dependence of the PISN mass scale depends on several astrophysical parameters (discussed above) and cannot be characterized uniquely. We show the inferred redshift of the value assuming a fixed PISN mass scale indicated by zPISN and the true redshift of the sources in the mock in Fig. 6 along with the uncertainty due to the inferred mass distribution for 2.5 years of LVK observations and 0.5 years of Cosmic Explorer for the case with t min d = 0.5 Gyr, κ = 1, γ = −0.44, α = 1.5. This scenario is chosen as it has the minimum shift in the PISN mass scale in comparison to the other cases. So, it will describe the minimum error for the scenarios considered in this analysis. The deviation is much larger for the third-generation detectors such as Cosmic Explorer and Einstein Telescope than for the current generation GW detectors such as LVK. The systematic error in the inferred redshift is going to vary for different choices of the astrophysical parameters t min d , κ, γ, α considered in this analysis. For the cases with large variation in the MPISN (see Fig. 3), the corresponding systematic error in the inferred redshift zPISN obtained using assuming a fixed of MPISN will be larger than the case shown in Fig. 6. The difference between the inferred redshift and true redshift is shown in Fig. 6 bottom panel for four different cases. The error bars on the redshift are the 68% C.I. obtained from the inferred redshift by combining the mass posteriors. This indicates that with the increase in more events, the statistical error in the mass inference will be subdominant than the systematic error due to the unknown redshift dependence of the PISN mass scale. We find that there is going to be about 10% variation for this model parameters. The variation can go as large as 25 − 30% for different values of the parameters t min d , κ, γ, α considered in this analysis (see Fig. 2). Moreover, this systematic uncertainty is only a minimum uncertainty. Incorrect astrophysical modeling, variation of metallicity with galaxies, multiple delay time distributions, and inaccurate prescriptions of stellar winds can increase this even further.

IMPACT ON GW SOURCE POPULATION INFERENCE AND COSMOLOGY
From the above discussion, it is evident that the redshift inference from the PISN mass scale is not accurate and cannot be well predicted due to two major sources of uncertainties (i) redshift dependence of the metallicity and (ii) redshift dependence of the delay time distribution. Even in the scenario when the PISN mass scale can be correctly modeled 2 , then also the observed mass population will not be redshift independent. This will have a significant impact on the inference of the GW source population as well as on the inference of cosmological results. Impact on the GW source population : The redshift dependence of the GW mass distribution will have a significant impact on understanding the true mass distribution of the sources. Our study shows that the position of the PISN mass cutoff and the mass range over which the distribution can be smeared both are going to depend on redshift through the delay time distribution and metallicity. As a result, it is important to include the redshift dependence of the mass distribution of the GW sources in the analysis. Moreover, the redshift dependence of the mass distribution of the BBHs and also their merger rates are going to be correlated due to their dependence on the delay time distribution. So, a joint analysis of the GW mass distribution and the merger rates will be appropriate to infer the PISN mass scale, metallicity dependence, and the delay time distribution. In a recent work (Karathanasis et al. 2022), we made a joint estimation of these quantities from the GWTC-3 data of the LVK collaboration. An independent way to infer the PISN mass distribution is using the cross-correlation technique after marginalizing the GW bias parameters (Mukherjee & Wandelt 2018;Mukherjee et al. , 2021bDiaz & Mukherjee 2021). A proper inference of the mass distribution of the BBHs will also be useful for the analysis of the stochastic GW background Callister et al. 2020b;Mukherjee & Silk 2021;Mukherjee et al. 2021a) and lensing event rates (Oguri 2018;Mukherjee et al. 2021e). The amplitude of the stochastic GW background and also the lensing event rates strongly depend on the mass distribution of the BBHs. Inclusion of the redshift dependence of the PISN mass distribution explored in this analysis will have a vital role in correctly estimating these signals from the data. Moreover, constraints on the lensing event rates using the stochastic GW background (Mukherjee et al. 2021c;Buscicchio et al. 2020) also need to use a mass distribution including the redshift dependence of the PISN mass-scale and the maximum mass.
Impact on the cosmological results : The redshift dependence of the PISN mass cutoff leads to a systematic bias in the inference of the true cosmological redshift. As shown in the previous section, the bias in the redshift estimation is mainly more severe at the high redshift than at the low redshift due to the dependence on metallicity. As shown in Fig. 6, the error in the redshift estimation can be about 10%. This error can be as large as 30% for the values of the parameters (t min d , κ, γ, α) considered in this analysis. This error is going to be subdominant from the statistical error for the current LVK data GWTC-3 (Abbott et al. 2021d,a,b,c), but in the future to achieve a few percent measures of the cosmological parameters such as the dark energy equation of state You et al. 2021), the systematic uncertainty from the PISN mass cutoff will be a bottleneck. In Fig. 7 we show the impact of the incorrect redshift inference on the luminosity distance (d ) and redshift plane for LVK and CE, assuming the value of MPISN as the one measured from low redshift (z = 0.2). We also show the true injected value in the yellow marker. To compare this with different expansion history models, we vary a single parameter w0 from -0.5 to -1.48, keeping all other cosmological parameters at a fixed value. This shows that for the inferred wrong redshifts using the PISN mass scale, there is going to be a significant bias in the value of inferred dark energy equation of state parameter towards a higher value (w0 > −1). This value is about a factor of two away from the true value. Similar to w0, the systematic effect will also influence the measurement of the Hubble constant and other cosmological parameters that appears in the Hubble parameter H(z). So, without carefully mitigating the systematic effect due to the PISN mass scale, one cannot use it for accurate inference of the cosmological parameters. The errors shown in Fig. 6 are for a scenario with a very simplistic model of the metallicity dependence and delay time distribution. The ignorance of the true values of these quantities leads to a systematic overestimate of the redshift of the GW sources. To reach a few percent level measurements of the cosmological parameters using the PISN mass cutoff, we need to mitigate this systematic error, otherwise, it will lead to a systematic bias in the inference of the cosmological parameters. Careful calibration of the mass distribution of BBH is required in order to mitigate possible systematic effects.
Similar to the background cosmological parameters, the use of the PISN mass cutoff was also thought to be useful to test the general theory of relativity (GR) from GW propagation by inferring the frictional term (Belgacem et al. 2018a,b;Ezquiaga 2021;Mancarella et al. 2021) and also for measuring extra dimensions (Hernandez 2021). The measurement of the deviation from the GR depends completely on the difference in the luminosity distance inferred from GW and the luminosity distance for photons at that redshift. Any systematic error in the redshift inference for the GW sources will lead to a bias in the value of frictional terms. For most of the non-GR theories, the deviation from GR is expected of the order of a few percent. As a result, to make a robust measurement of any deviation from GR, it will be essential to correctly include the systematic error that can be induced from the The difference between the inferred redshift z PISN and the true redshift ztrue is shown for LVK with 2.5 years and 10 years and for one CE with 0.5 years and 5 years of observation time.
redshift dependence of the PISN mass cutoff. Our analysis shows that there can be a 10-30% systematic error in PISN mass scale over the redshift up to z = 7. As a result, a few percent level measurement is not possible unless this effect can be properly mitigated. Several recent analyses (Ezquiaga 2021;Mancarella et al. 2021;Hernandez 2021) assume a fixed PISN mass scale and hence these measurements are not reliable and subject to systematic errors. The inference of the cosmological parameters jointly with the frictional term from GW propagation can be made in a robust way using the cross-

CONCLUSION AND THE WAY FORWARD
In this work, we explore for the first time the redshift dependence of the PISN mass cutoff and its implication on the estimation of the source redshift of BBHs. We show that even if the PISN mass cutoff and its dependence on the stellar metallicity can be well modeled theoretically, the redshift evolution of the metallicity and the delay time distribution of the BBH mergers can lead to a mixing between BHs formed across a large cosmic epoch. This will lead to redshift dependence of the PISN mass cutoff. The unknown form of the probability distribution of the delay time will lead to an unknown amount of systematic shift in the PISN mass cutoff to higher values with a change in the redshift. The mass scales over which the BH distribution will be smeared to zero also depend on the metallicity and the delay time distribution. Ignoring the redshift dependence for any cosmological and GW population analysis can lead to a systematic bias.
We show that the use of a fixed PISN mass cutoff to infer the source redshift of the BBHs can lead to a biased inference. Even for one of the best cases (minimum shift in the PISN mass cutoff) and assuming that one can model the metallicity dependence of the PISN mass cutoff theoretically, we find that there is going to be around 10 − 30% systematic error in the inference of the redshift of the sources (towards high values). Though this error is sub-dominant to the statistical error in the mass measurement of the current LVK data (Abbott et al. 2021d,a,b,c) and the cosmological results inferred from GWTC-3 (Abbott et al. 2021c) are unlikely to be influenced by this, but it is going to be a major source of uncertainty in the future when statistical error will be subdominant. This is also in agreement with the recent studies from GWTC-2 and GWTC-3 (Fishbach et al. 2021;Abbott et al. 2021a). We show that for a network of four GW detectors such as LIGO-Hanford, LIGO-Livingston, Virgo, and KAGRA and also in the future from Cosmic Explorer, the systematic error will lead to a systematic bias in the inferred value of redshift at the level of 10-30%. As a result, to reach a percent level measurement of the cosmological parameter using the PISN mass cutoff, it is essential to mitigate this systematic error. Similarly, a wrong inference of the redshift will have an impact on the measurement of any non-GR signatures due to GW propagation.
The redshift dependence of the PISN mass cutoff modeled in this paper will also be useful for understanding the source frame mass distribution of the BBHs. The inclusion of a redshift dependence of the PISN mass cutoff and exploring its dependence on the stellar metallicity and delay time distribution will able to shed light on the mass distribution of the BHs in the Universe. A systematic shift in the PISN mass cutoff with redshift and its dependence on the delay time will be able to help us in understanding the formation channel of the BBHs. An independent way to infer the redshift dependence of the PISN mass scale is by inferring the clustering redshift using the cross-correlation technique and marginalizing over the GW bias parameters (Oguri 2018;Mukherjee & Wandelt 2018;Mukherjee et al. 2021b;Scelfo et al. 2020;Diaz & Mukherjee 2021).
In summary, this work shows for the first time the redshift independence of the PISN mass cutoff is not possible in reality. Even when there is no theoretical modeling uncertainty, astrophysical uncertainties will play a significant role in inducing the redshift dependence of the PISN mass cutoff. In addition, there can be also secondary sources, and additional theoretical uncertainties that even increase the budget of the systematic error. As a result, the use of the PISN mass cutoff for inferring redshift to the source, and use that for cosmological analysis can cause a significant systematic bias. Though currently, the statistical error is larger than the systematic error, in the future, an accurate measurement of the cosmological parameters will not be possible unless this can be mitigated by other techniques. In future work, we will develop techniques to also include the redshift dependence of the PISN mass cutoff and its joint estimation with the cosmological parameters.