Modified gravitational wave propagation and the binary neutron star mass function

Modified gravitational wave (GW) propagation is a generic phenomenon in modified gravity. It affects the reconstruction of the redshift of coalescing binaries from the luminosity distance measured by GW detectors, and therefore the reconstruction of the actual masses of the component compact stars from the observed (`detector-frame') masses. We show that, thanks to the narrowness of the mass distribution of binary neutron stars, this effect can provide a clear signature of modified gravity, particularly for the redshifts explored by third generation GW detectors such as Einstein Telescope and Cosmic Explorer.


Introduction
In recent years, modified gravitational wave (GW) propagation has come to attention as one of the most promising ways of testing deviations from General Relativity (GR) on cosmological scales. The effect is encoded in the propagation equation of GWs across cosmological distances which, in modified gravity theories, can take the form [1][2][3][4][5][6][7][8] whereh A (η; k) is the Fourier transform of the GW perturbation, h = ∂h/∂η where η is conformal time, a(η) is the scale factor, H = a /a, and A = +, × labels the two polarizations. The difference with respect to GR is given by a nonvanishing function δ(η). Several other modifications with respect to GR are possible in the propagation equation of GWs. The most immediate options are a deviation of the speed of GWs from the speed of light, or a graviton mass. Both would rather modify the c 2 k 2 term in eq. (1), but are now very significantly constrained: a deviation of the speed of GWs from the speed of light is excluded at the level |c gw − c|/c < O(10 −15 ) by the observation of GW170817 and its electromagnetic counterpart [9] (and a large class of modifications of GR have been ruled out by this limit [10][11][12][13]), while limits on the graviton mass are in the range O(10 −32 − 10 −22 ) eV, depending on the probes used [14]; several other modifications, in general related to rather specific classes of modified gravity theories, have been tested or proposed, such as extra polarizations [15], Lorentzviolating dispersion relations [16], parity-violating effects [17], or scale dependent modifications of the speed of GWs [18].
Modified GW propagation, in the form described by eq. (1), was first found in some explicit scalar-tensor theories of the Horndeski class [1][2][3][4] (see also [19] for a discussion within the effective field theory approach to dark energy) and, in refs. [5,7], in non-local infrared modifications of gravity, i.e. in theories where the underlying classical action is still GR, but non-local terms, relevant in the infrared, are assumed to be generated by non-perturbative effects in the quantum effective action [20] (see [21] for recent review). However, it has been understood that the phenomenon is completely general and appears in all best studied modified gravity theories [8]. It also appears, in a different form not described by eq. (1), in theories with extra dimensions, where it is rather due to the loss of gravitons to the bulk [22,23].
The modified friction term in eq. (1) changes the evolution of the GW amplitude in the propagation across cosmological distances. Since, in GR, the amplitude of a coalescing binary is proportional to 1/d L , where d L is the luminosity distance, this introduces a bias in the luminosity distance inferred from GW observations. In particular, if δ(η) < 0, the damping term is stronger and, after propagation from the source to the detector, the GW has a smaller amplitude. If interpreted within GR, it would therefore appear to come from a distance larger than its actual distance (and vice versa for δ(η) > 0). It is then useful to introduce a distinction between the standard luminosity distance, that, in this context, is called the 'electromagnetic luminosity distance' and denoted by d em L , and the luminosity distance extracted from the observation of the GWs of a compact binary coalescence, that is called the 'GW luminosity distance' [5] and denoted by d gw L . The two quantities are related by [5,7] d gw where δ(z) ≡ δ[η(z)]. A useful parametrization of this effect, which catches the redshift dependence predicted by almost all explicit models in terms of just two parameters (Ξ 0 , n), is ob-tained writing [7], which interpolates between d gw L /d em L = 1 as z → 0 and an asymptotic value Ξ 0 at large z, with a power-law behavior in a = 1/(1 + z) fixed by n. GR is recovered when Ξ 0 = 1 (for all n). The study of explicit modified gravity models shows that Ξ 0 can be significantly different from 1. In particular, in nonlocal gravity it can be as large as 1.80 [21,24], corresponding to a 80% deviation from GR, despite the fact that this model complies with existing observational bounds, that force deviations from GR and from ΛCDM in the background evolution and in the scalar perturbation sector to be at most of a few percent [25,26]. Thus, the newly opened window of GWs could give us the best opportunities for testing modified gravity and dark energy.
Contrary to quantities such as the speed of GWs or the graviton mass, the limits on the parameter Ξ 0 (the main parameter that describes modified GW propagation; the power n in eq. (3) only determines the precise form of the interpolation between the asymptotic values) are still quite broad. Using the binary neutron star (BNS) GW170817, with the redshift determined from the electromagnetic counterpart, only gives bounds of order Ξ 0 < ∼ 14 (68% c.l.) [7] (see also [4,27]). This is because the redshift of GW170817 is very small, z 0.01, and d gw L (z)/d em L (z) goes to one as z → 0, for all Ξ 0 . A more significant limit, has been obtained in [28], using binary black hole (BBH) coalescences without electromagnetic counterpart ('dark sirens') from the O1, O2 and O3a runs of the LIGO/Virgo Collaboration (LVC) and correlating them with the GLADE galaxy catalog [29]. An even more stringent measurement is obtained under the tentative identification of the flare ZTF19abanrhr as the electromagnetic counterpart of the BBH coalescence GW190521, in which case one gets Ξ 0 = 1.8 +0.9 −0.6 [28] (see also [30]). However, this identification currently is not secure. A limit on modified GW propagation (using a different parametrization) has been obtained in [31] using the BBH mass function, following an idea originally proposed in [32] to infer H 0 , and a recent re-analysis in [33], using again the BBH mass function, gives while the corresponding limit at 90% c.l. is Even with the study based on the BBH mass function, which currently gives the most stringent bounds on Ξ 0 , current data are not constraining enough to obtain a limit on n, with the posterior reflecting basically the prior used [33].
Since the effect of modified GW propagation increases with redshift (at least until the ratio in eq. (3) saturates to its large z limit d gw L /d em L Ξ 0 ), third generation (3G) ground based GW detectors such the Einstein Telescope (ET) [34,35] and Cosmic Explorer (CE) [36], or the space interferometer LISA [37], are particularly well suited to study it, and several forecasts have been made on the accuracy that future observations can reach on Ξ 0 , using different techniques [8,[38][39][40][41][42][43][44][45][46]. Observe that the function δ(η) in eq. (1) only affects the amplitude of the GW signal. Other effects, such as modified dispersion relations, also affect the post-Newtonian coefficients of the phase (see [47] for the most recent bounds using LIGO/Virgo data). In that case, eventually, a joint analysis of the Hubble parameter H 0 , of modified GW propagation and of modified dispersion relations might be necessary [40]. 1 The same logic that has been used in [31,33] to obtain bounds on modified GW propagation from the BBH mass function can be applied to the BNS mass function, with the advantage that the latter is narrow and is not expected to evolve significantly with redshift. The idea of using the BNS mass function for extracting cosmological informations was proposed in [48][49][50] in the context of the determination of H 0 within ΛCDM. In this paper we will discuss its application to modified GW propagation. We will see that the method based on the BNS mass function can become particularly powerful when applied to modified GW propagation at 3G detectors, thanks to the fact that modified GW propagation increases with distance, and ET and CE can detect BNS up to large redshifts, z 2 − 3 for ET, and even z ∼ 10 for CE [51]. 2

Modified GW propagation and mass reconstruction
The starting point of our analysis is the fact that GW detectors measure the GW luminosity distance of the source, d gw L , which is different from the actual electromagnetic luminosity distance d em L if, in Nature, Ξ 0 1. The redshift z GR of the source inferred from the measured d gw L assuming GR and ΛCDM (with a given value of H 0 and Ω M , that we keep for definiteness the same in GR and in the modified gravity theory under consideration), would therefore differ from the true value z true . The effect is shown in Fig. 1, for a sample of values of Ξ 0 consistent with eq. (6). We see that the effect can become very significant at large redshifts.
In turn, this affects the reconstruction of the actual masses m i (i = 1, 2) of the component stars ('source-frame' masses, as they are called in this context), from the 'detector-frame' masses m (det),i ≡ (1 + z)m i , that are the quantities directly obtained from the GW observations. If Nature is described by a modified gravity theory with Ξ 0 1, the true values of the source-frame masses, m true,i , are related to the values of the source-frame masses that would be inferred in GR, m GR,i , by where m GR,i = m (det),i /(1 + z GR ). The same multiplicative bias factor will appear in any other combination with dimensions of mass of the source-frame masses of the component stars, such as the total source-frame mass m tot = m 1 + m 2 , or the sourceframe chirp mass M c = (m 1 m 2 ) 3/5 /m 1/5 tot . The upper panel of Fig. 2 shows the ratio M GR /M true , as a function of z true (while the lower panel shows M true /M GR , as a function of z GR ), for any such mass scale. We see that, at the redshifts accessible to ET and CE, and for values of Ξ 0 consistent with current limits, the effect can be very large. For instance, setting Ξ 0 = 1.8, for a NS with m true = 1.35M at z true = 1 (that, with this value of Ξ 0 , corresponds to z GR 1.45), the mass incorrectly inferred from GR would be m GR 1.10M ; at z true = 2 (z GR 3.10) for the same system in GR one would infer m GR 0.99M ; and, for a BNS with the same mass at z true = 5 (z GR 8.20), which could still be accessible to CE, one would find m GR 0.88M . Furthermore, exactly the same factor affects the two component stars (which is not the case in general for astrophysical effects), so a BNS with (1.35 + 1.35)M would appear as a (1.10 + 1.10)M system for z true = 1, as a (0.99 + 0.99)M system for z true = 2, and as a (0.88 + 0.88)M system at z true = 5. Compared to the narrowness of the neutron star (NS) mass distribution, this is a huge effect. The mass of the BNSs found with electromagnetic observations can be described by a Gaussian distribution with mean 1.33M and standard deviation 0.09M [52] (which, assuming that the distribution of the two masses are independent, corresponds to a Gaussian distribution for the total mass with mean 2.66M and standard deviation 0.13M ), or by a flat distribution between a minimum and a maximum mass, with a similar width. Somewhat broader limits are obtained from an analysis using only the NSs in BNS or in BH-NS systems detected by GW observations [53], although this sample, of six NSs, is very limited.
These estimates show that even a single BNS at large z would have a significant constraining power on Ξ 0 . Still, if one would find just a single system that, interpreted within GR, corresponds to, say, a (1.0 + 1.0)M binary at z GR 3.1, as in one of the examples above, one would remain in doubt on whether this is a binary made of exotic compact objects, such as primordial black holes, or a signal of modified GW propagation. The power of the method, however, is that the same effect will affect all BNS systems, by a factor that depends only on z. If Nature is described by a modified gravity theory with a large deviation from GR such as, say, Ξ 0 = 1.8, as in the examples above, at large redshifts ET and CE will not find a single BNS whose component masses, interpreting the data within GR, will be near the typical value of 1.35M . When interpreted within GR, all BNS with z true = 1 would appear to have component masses around 1.10M ; all BNS at z true = 2 would appear to have masses around 0.99M , and so on. The detection rate of BNS at ET and CE will be impressive, of order of 7×10 4 events per year already for a single detector such as ET [38,54,55] and, among these, within a GR interpretation, there would not be a single 'normal' neutron star at large z, but rather a plethora of objects with puzzling masses. The situation is illustrated in Fig. 3, where m GR tot denotes the total (source-frame) mass of the BNS inferred in GR, for different values of Ξ 0 . Here we have assumed that the distribution of the source-frame total mass of the binary is a Gaussian, with mean 2.66M and standard deviation 0.13M . In the absence of astrophysical evolutionary effects, for which, currently, there is little observational information, but which are not expected by any means to give ef- fects comparable to those shown in the figure (see, e.g., Fig. 4 of ref. [56]), the distribution would not change with redshift (black dotted line). In the presence of modified GW propagation, with the values of Ξ 0 shown in the figure, that represent deviations from GR large but still consistent with current limits, the masses wrongly inferred using GR are narrowly distributed around completely different mean values. Note that the apparent skewness of the Gaussians at large z is just a graphical effect in this plot (selection effects could, however, introduce some actual skewness, since the low mass end will be less detected at high redshifts).
Actually, it is convenient to use the chirp mass, rather than the total mass, because, in GW observations, the chirp mass is measured much more accurately than the individual component masses or the total mass. The corresponding apparent evolution in redshift is shown in Fig. 4. Actually, the distribution of total masses, and also of chirp masses, can be fitted both by a Gaussian distribution of by a flat distribution between minimum and maximum values. In Fig. 4 we have used, for illustration, a flat distribution between M c,min = 1.10M and M c,max = 1.25M , that encompasses the chirp masses of all BNS that merge within a Hubble time, reported in Table 1 of [52].
Finally, another important signature of modified GW propagation will be given by how the BNS population is distributed in redshift (i.e., the absolute normalization of the distributions, that in Fig. 3 and 4 have been normalized to unity). Even if our prior information on the BNS merger rate is not as stringent as on the BNS mass function, still we expect that the rate will be where C 0 (z p , α z , β z ) ≡ 1 + (1 + z p ) −α z −β z is a normalization constant that ensures R(0) = R 0 , and z p is the peak of the star formation rate, which is known to be in the range z p (2 − 3). In a modified gravity theory, the difference between z GR and z true will lead to a bias in the reconstruction of R(z). 3 For instance we have seen that, if Nature is described by a modified gravity theory with our reference value Ξ 0 = 1.8, and we rather use GR to interpret the data, a BNS with z true = 2 would be wrongly interpreted as having a redshift z GR 3.10, and z true = 3 corresponds to z GR 4.79. The peak of the BNS merger distribution would then appear to be at redshifts larger than the peak of the star formation rate, leading to another puzzling result of the GR interpretation (that, for Ξ 0 > 1, could not be explained in terms of delay between formation and merger, since in this case one would find that the peak of the merger rate took place before the peak of the star formation rate). A joint Bayesian inference on the BNS mass function and on the BNS rate parameters would therefore further strengthen the power of the method.

Sources of errors
The above discussion is still idealized, because it neglected the errors on the measurements. The relative accuracy on the detector-frame chirp mass M c = (1 + z)M c is of order ∆M c /M c ∼ 1/N c , where N c is the number of inspiral cycles of the signal in the detector bandwidth, see, e.g. eq. (7.187) of [60]. For a lower cutoff of the detector near 3 Hz, as in the design of ET, and the chirp mass of a BNS, we have N c 10 5 (using eq. (4.23) of [60]). The error on the detector-frame chirp mass is therefore negligible. More important is the error on the redshift due to the observational error on d The function f (z; Ξ 0 ) goes from zero at z = 0 to one at large z, with only mild dependence on Ξ 0 . The value of ∆d gw L /d gw L as a function of redshift can be obtained using the fitting formulas provided in [38] [eq. (2.13) for ET, and eq. (2.20) for a network ET+CE+CE], which were obtained from a mock source catalog of BNS detections, averaging over detector orientation, source inclination, and BNS mass distribution. In Fig. 5 we show the result for ∆M c /M c at ET (upper panel) and at a network ET+CE+CE (lower panel). We see that, on average, the relative error on the source-frame chirp mass, induced by the observational error on d gw L , is below 6% up to z < ∼ 9 for a network ET+CE+CE (where ET contributes to BNS detections only up to z 3). Similarly, we find that it is below (5 − 6)% up to z < ∼ 3 for ET alone. This is smaller than the intrinsic relative width of the BNS mass distribution, ∆m/m ∼ 0.1 obtained from electromagnetic observations of BNS, and therefore also of the corresponding distribution of chirp masses. So, the accuracy of the method appears to be limited more by the intrinsic width of the BNS mass distribution, than by observational errors on the reconstruction of the redshift.
Other sources of error would require more complex dedicated studies, that are beyond the scope of this paper. One is the error due to lensing from large scale structures along the line of sight. On linear scales, inhomogeneities induce a relative error ∆d L /d L < ∼ 1% for all redshifts z < 5 [61] (see also Fig. 12 of [35]). Therefore, these are smaller than the error on the measurement of the luminosity distance in ET. The treatment of non-linear scales is, however, more complex and has been recently discussed in [62]. In this case, using a simplified model for modified GW propagation, corresponding to setting δ(z) equal to a constant δ 0 in eq. (2), 4 and using the notation ν = −2δ 0 , ref. [62] finds that, at ET, 350 BNS events with counterpart would be needed to measure ν at the 1% level. It would be very interesting to extend the study in [62] to the computation of the effect of lensing from clustered structures on the reconstructed BNS mass function, using furthermore the full expression (2) for modified GW propagation.
Another important point concerns the evolution with redshift of the BNS mass function, due to astrophysical evolutionary effects. For values of Ξ 0 such as our reference value Ξ 0 = 1.8, that represents a large deviations from GR, we have seen in Figs. 3 or 4, and in the discussion around them, that the effect of modified GW propagation is very large compared to anything that could be expected from evolutionary effects. We certainly do not expect that at, say, redshift z = 2, neutron stars have a mass m 1.0M , and they only reach the typical observed values m 1.33M in the local Universe because of evolutionary effects; in contrast, as we have seen, m 1.0M is the value that would be erroneously reconstructed using GR, if Ξ 0 = 1.8 and z true = 2. So, at this level, evolutionary effects cannot mimic modified gravity. However, to detect finer deviations from GR, corresponding to values of Ξ 0 closer to one, eventually also evolutionary effects in BNS will have to be taken into account. By the time that ET and CE will be operational, more information on the evolution of the BNS masses with redshift might have been obtained from electromagnetic observations, expanding the currently very limited sample of 17 BNS, used in [52]. Furthermore, as we already remarked, a handle to discriminate the effect of modified GW propagation from evolutionary effects is that the former acts exactly in the same way on the reconstruction of the two star masses, and therefore does not affect the inferred mass ratio, while evolutionary effects in general modify the mass ratio. Eventually, the best strategy will be to perform a joint inference of the cosmological parameters and of the parameters describing the astrophysical population, along the lines discussed in [31][32][33][63][64][65], including all this prior information.
A full Bayesian analysis on mock data for 3G detectors, including selection effects and the current understanding of observational errors, which is necessary to reliably quantify the accuracy that can be obtained on Ξ 0 , is under development and will be presented in a separate paper.

Conclusions
For BNSs, at the large redshifts that will be probed by thirdgeneration detectors such as Einstein Telescope and Cosmic Explorer, modified GW propagation could leave a very characteristic imprint on the mass distribution (and, to some extent, also on the redshift distribution) of the observed BNS. In modified gravity, the size of the effect is controlled by the parameter Ξ 0 introduced in eq. (3). For values as large as Ξ 0 1.8, that are consistent with current limits and are on the upper range of the predictions from an explicit and viable model of modified gravity [20,21,24], the effect on the reconstructions of the BNS mass function and of the BNS merger rate is quite striking, and would provide a clear and unambiguous signature of modifications of General Relativity on cosmological scales. For values of Ξ 0 closer to the GR value Ξ 0 = 1, disentangling the effect of modified gravity from astrophysical and cosmological effects (such as evolutionary effect in the BNS mass function or lensing by non-linear structures) will eventually become more challenging and will require further studies.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.