Impacts of gravitational-wave standard siren observations from Einstein Telescope and Cosmic Explorer on weighing neutrinos in interacting dark energy models

The multi-messenger gravitational-wave (GW) observation for binary neutron star merger events could provide a rather useful tool to explore the evolution of the universe. In particular, for the third-generation GW detectors, i.e., the Einstein Telescope (ET) and the Cosmic Explorer (CE), proposed to be built in Europe and the U.S., respectively, lots of GW standard sirens with known redshifts could be obtained, which would exert great impacts on the cosmological parameter estimation. The total neutrino mass could be measured by cosmological observations, but such a measurement is model-dependent and currently only gives an upper limit. In this work, we wish to investigate whether the GW standard sirens observed by ET and CE could help improve the constraint on the neutrino mass, in particular in the interacting dark energy (IDE) models. We find that the GW standard siren observations from ET and CE can only slightly improve the constraint on the neutrino mass in the IDE models, compared to the current limit. The improvements in the IDE models are weaker than in the standard cosmological model. Although the limit on neutrino mass can only be slightly updated, the constraints on other cosmological parameters can be significantly improved by using the GW observations.


I. INTRODUCTION
In the recent two decades, the study of cosmology has entered the era of precision cosmology. A standard model of cosmology has been established, usually called the Λ cold dark matter (ΛCDM) model. The measurements of cosmic microwave background (CMB) anisotropies from the Planck satellite mission have constrained the six primary parameters of the ΛCDM model with unprecedented precision. However, with the measurement precisions of the cosmological parameters improved, some puzzling issues appeared. For example, the inferred values of the Hubble constant from the Planck observation of the CMB anisotropies (based on the ΛCDM model) [1] and from the Cepheid-supernova distance ladder measurement [2] are inconsistent, with the tension between them more than 4σ significance [2]. Namely, there is an inconsistency of measurements between the early and late universe, which is the so-called "Hubble tension" problem. The Hubble tension recently has been widely discussed in the literature (see, e.g., Refs. ). Furthermore, theoretically, for the ΛCDM model, the cosmological constant Λ, which is equivalent to the density of vacuum energy, has always been suffering from serious theoretical challenges, such as the "fine-tuning" and "cosmic coincidence" problems [32,33]. Thus, it is hard to say that the ΛCDM model with only six base parameters is the eventual model of cosmology. All of these facts actually imply that the ΛCDM model needs to be further extended and some extra parameters concerning new physics need to be introduced into the new models. Of course, some novel cosmological probes should also be further developed.
To extend the ΛCDM cosmology, the primary idea is to consider the dynamical dark energy with the darkenergy density no longer a constant. In this class of models, the simplest one is the model with a dark energy having a constant equation-of-state (EoS) parameter, w = p de /ρ de = constant, which is usually called the wCDM model. For some popular dark energy models, see, e.g., Refs. [34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51]. There is also a class of models known as the interacting dark energy (IDE) models in which some direct, non-gravitational interaction between dark energy and dark matter is considered. The interaction between dark sectors could help resolve (or alleviate) the coincidence problem of dark energy, and also can help alleviate the Hubble tension. The IDE models have been widely studied and deeply explored till now (see, e.g., Refs. [11,13,).
Currently, the mainstream cosmological probes mainly include, e.g., the CMB anisotropies, baryon acoustic oscillations (BAO), type Ia supernovae (SN), direct determination of the Hubble constant (H 0 ), weak gravitational lensing, redshift space distortions, and clusters of galaxies. The combinations of these cosmological data based on the electromagnetic (EM) observations have provided precise measurements for the base cosmological parameters. But for some extended parameters beyond the standard ΛCDM model, e.g., the EoS parameter of dark energy, the tensor-to-scalar ratio, and the total mass of neutrinos, they still cannot be precisely measured. Therefore, on one hand, the EM observations should be further greatly developed, and on the other hand, some novel cosmological probes are also needed to be developed in the future. In the next few decades, the gravitationalwave (GW) standard siren observation is one of the most promising cosmological probes.
The detection of the GW event GW170817 [98] from the binary neutron star (BNS) merger initiated the multimessenger astronomy era. Because in this event, not only GWs but also the EM signals in various bands were detected for the same transient source [99,100]. From the analysis of the waveform of GW, one can obtain the absolute luminosity distance to the source. Furthermore, the redshift of the source can also be determined by identifying the EM counterpart of the GW source. With the known distance and redshift of a celestial source, a distance-redshift relation can be established, which can be used to explore the expansion history of the universe [101]. Such a tool of exploring the evolution of the universe provided by GWs is called "standard sirens" (note here that the case having EM counterparts is usually referred to as bright sirens, to be differentiated with the case of dark sirens without EM counterparts) [102].
The main advantage of the GW standard siren method is that the absolute luminosity distances can be measured. This is obviously superior to the SN observation, in that the latter can only measure the ratio of luminosity distances at different redshifts. In addition, the GW observation can observe much higher redshift events, compared to the SN observation.
It is indisputable that the GW standard siren will be developed into a powerful cosmological probe in the future. The third-generation (3G) ground-based GW detectors, such as the Einstein Telescope (ET) [103,104] in Europe and the Cosmic Explorer (CE) [105,106] in the United States, have been proposed. In the 2030s, ET will be brought into operation. CE will start its observation in the mid-2030s. The 3G ground-based GW detectors have a much wider detection-frequency range and a much better detection sensitivity, which can observe much more BNS events at much deeper redshifts. Recently, the GW standard sirens have been widely discussed in the literature  (see Ref. [135] for a recent review). It is found that the GW standard siren observations from ET and CE would play an important role in the cosmological parameter estimation [136][137][138][139][140][141].
In cosmology, neutrinos play a crucial role in helping shape the large-scale structure and the expansion history of the universe. The phenomenon of neutrino oscillation indicates that neutrinos have masses and there are mass splittings between different neutrino species. However, it is extremely difficult to measure the absolute masses of neutrinos. Neutrino oscillation experiments cannot measure the absolute neutrino masses, but can only give the squared mass differences between the different mass eigenstates of neutrinos. The solar and reactor experiments give ∆m 2 21 7.5 × 10 −5 eV 2 and the atmospheric and accelerator beam experiments give |∆m 2 31 | 2.5 × 10 −3 eV 2 [142,143]. Thus, there are two possible mass hierarchies of the neutrino mass spectrum, namely, the normal hierarchy (NH) with m 1 < m 2 m 3 and the inverted hierarchy (IH) with m 3 m 1 < m 2 . In addition, in some cases one also considers the cosmological models of neglecting the neutrino mass splittings, namely m 1 = m 2 = m 3 , which is usually called the degenerate hierarchy (DH).
Although the neutrino masses can hardly be measured by particle physics experiments, they can be effectively constrained by the cosmological observations. This is because massive neutrinos can exert some impacts on the evolution of the universe. Using the current cosmological observations, an upper limit on the total neutrino mass m ν can be obtained. So far, the most stringent limit on the total neutrino mass comes from the Planck 2018 CMB observation, and the combination CMB+BAO+SN gives the 95% CL upper limit m ν < 0.12 eV, for the DH case in the ΛCDM model. See e.g. Refs.  for studies on neutrino mass in cosmology.
In a recent forecast [138], it was shown that the standard sirens observed by the ET can be used to improve the constraints on the total neutrino mass in the ΛCDM model. Using 1000 GW standard siren data points of the BNS merger events, it is found that the upper limits on m ν can be tightened by about 10% [138]. However, weighing neutrinos in cosmology depends on the cosmological model considered, and thus one would be curious about whether the role the GW data play in helping measure the neutrino mass will change if an extension to the ΛCDM model is considered. In this work, we consider the IDE models, and we wish to see what will happen on measuring neutrino mass when the IDE models are considered.
In an IDE model, the energy conservation equations for dark energy and CDM satisfẏ where Q is the energy transfer rate, ρ de and ρ c represent the energy densities of dark energy and CDM, respectively, H is the Hubble parameter, and a dot represents the derivative with respect to the cosmic time t. In this work, we consider the interaction form of Q = βHρ c , where β is a dimensionless coupling parameter. Here, β > 0 and β < 0 means CDM decaying into dark energy and dark energy decaying into CDM, respectively. In this work, we consider the IDE versions of the ΛCDM and wCDM models, which are called the IΛCDM and IwCDM models. We will discuss the cosmological parameter estimation in the IΛCDM+ m ν and IwCDM+ m ν models. Moreover, we will consider the three neutrino mass hierarchy cases, i.e., the NH, IH, DH cases. To avoid the perturbation divergence problem in the IDE models, in this work we employ the extended parameterized post-Friedmann (ePPF) framework [182,183] to calculate the perturbations of dark energy.
We simulate the GW standard siren data observed by ET and CE, and we use these simulated GW data to investigate how well they can be used to improve the constraints on the neutrino mass as well as other cosmological parameters on the basis of the current CMB+BAO+SN constraints.
The rest of this paper is organized as follows. In Sec. II A, we introduce the methods of simulating the GW standard siren data. In Sec. II B we describe the EM cosmological observations used in this work. In Sec. II C, we briefly describe the methods of constraining cosmological parameters. In Sec. III, we give the constraint results and make some relevant discussions. The conclusion is given in Sec. IV.

II. METHOD AND DATA
In this section, we first introduce the method of simulating the GW standard siren data from ET and CE. Then, we describe the current mainstream EM cosmological observations used in this work. Finally, we briefly introduce the method of constraining cosmological parameters.

A. Simulation of the GW standard sirens
The primary GW sources in the detection frequency band of the ground-based GW detectors are the mergers of BNS, binary stellar-mass black hole (BBH), and so on. The BNS mergers could produce rich EM signals [184] that can be detected by the EM observatories, thus enabling precise redshift measurements. Owing to the fact that there are no EM signals produced in the process of the BBH mergers, their redshifts could not be precisely measured through the detection of the EM counterparts. Hence, in this work, we only simulate the GW standard sirens from the BNS mergers.
Following Refs. [107,185], the redshift distribution of the BNS mergers takes the form where d C (z) is the comoving distance at the redshift z, and R(z) represents the redshift evolution of the burst rate, which takes the form [107,186,187] In Fig. 1, we show the redshift distribution of BNS mergers.
Considering the transverse-traceless gauge, the strain h(t) in the GW interferometers can be described by two independent polarization amplitudes, h + (t) and h × (t), where F + and F × are the antenna response functions, θ and φ describe the location of the GW source relative to the GW detector, and ψ is the polarization angle.
The antenna response functions of ET are [185] F (1) Since ET has three interferometers with 60 • inclined angles between each other, the other two response functions are F For CE, the antenna response functions are We consider the waveform in the inspiralling stage for a non-spinning BNS system. Here we adopt the restricted post-Newtonian (PN) approximation and calculate the waveform to the 3.5 PN order [188,189], where the Fourier amplitude A is given by where M c = (1 + z)M η 3/5 is the observed chirp mass, M = m 1 + m 2 is the total mass of the coalescing binary system with m 1 and m 2 being the component masses, η = m 1 m 2 /M 2 is the symmetric mass ratio, d L is the luminosity distance to the GW source, ι is the inclination angle between the binary's orbital angular momentum and the line of sight, t c is the coalescence time, and ψ c is the coalescence phase. The definitions of the functions Ψ and ϕ (2.0) can refer to Refs. [188,189].
The signal-to-noise ratio (SNR) for the network of N (i.e., N = 3 for ET and N = 1 for CE) independent interferometers can be calculated by where The inner product is defined as where f lower is the lower cutoff frequency (f lower = 1 Hz for ET and f lower = 5 Hz for CE), f upper = 2/(6 3/2 2πM obs ) is the frequency at the last stable orbit with M obs = (m 1 +m 2 )(1+z) [185], and S n (f ) is the oneside noise power spectral density (PSD). We adopt PSD of ET from Ref. [190] and that of CE from Ref. [191]. In this work, we choose the threshold of SNR to be 8 in our simulation.
For the 3G ground-based GW detectors, a few ×10 5 BNS merger events per year could be detected, but only about 0.1% of them may have γ-ray bursts toward us [192]. Recently, in Ref. [193] Chen et al. made a forecast and showed that 910 GW standard siren events could be observed by a 10-year observation of CE and Swift++. Therefore, in our forecast in the present work, for ET and CE, we simulate 1000 GW standard sirens from BNS mergers based on the 10-year observation.
We consider three measurement errors of d L , consisting of the instrumental error σ inst dL , the weak-lensing error σ lens dL , and the peculiar velocity error σ pv dL . Therefore, the total error of d L is We first use the Fisher information matrix to calculate σ inst dL . For a GW event, when using Fisher information matrix to estimate σ inst dL , we consider a 9 × 9 Fisher information matrix consisting of nine parameters of a GW source (d L , M c , η, θ, φ, ι, t c , ψ c , ψ), and thus the correlations between the nine parameters are considered in the analysis. For a network of N independent interferometers, the Fisher information matrix can be written as with h given by where θ i denotes nine parameters (d L , M c , η, θ, φ, ι, t c , ψ c , ψ) for a GW event. Then we have where F ij is the total Fisher information matrix for the network of N interferometers. Note that here σ inst dL = ∆θ 1 .
The error caused by weak lensing is adopted from Refs. [194,195], Here, we consider a delensing factor F delens . We use dedicated matter surveys along the line of sight of the GW event in order to estimate the lensing magnification distribution, which can remove part of the uncertainty due to weak lensing. This reduces the weak lensing uncertainty. The delensing factor is given by The error caused by the peculiar velocity of the GW source is given by [196] σ pv where H(z) is the Hubble parameter. v 2 is the peculiar velocity of the GW source and we roughly set v 2 = 500 km s −1 .
For each simulated GW source, the sky location, the binary inclination, the coalescence phase, and the polarization angle are randomly chosen in the ranges of cos θ ∈ [−1, 1], φ ∈ [0, 360 • ], ι ∈ [0, 20 • ], ψ c ∈ [0, 360 • ], and ψ ∈ [0, 360 • ]. The mass of an NS is randomly chosen in the range of [1,2] M . Without loss of generality, the merger time is chosen to t c = 0 in our analysis. Here we wish to note that the inclination angle should be randomly chosen in the range of cos ι ∈ [−1, 1] when simulating isotropic GW sources. However, in this work, we simulate GW events by detecting short γ-ray bursts (SGRBs) to determine sources' redshifts. Owing to the fact that SGRBs are strongly beamed [197], the detectable inclination angle is about ι ≤ 20 • [192,198]. Hence, in the present work, we set the inclination angle to be in the range of [0, 20 • ]. This is an ideal treatment, but for this work, since the number of simulated GW standard sirens is fixed, it has little effect on showing the impact of GW standard sirens on breaking cosmological parameter degeneracies and improving constraints on the cosmological parameters. In Fig. 2, we show the σ dL /d L scatter plot of the simulated standard siren data from ET (upper panel) and CE (lower panel). 1 We can see that: (i) the total errors 1 Here we wish to note that in the colorbars of Fig. 2, we set the maximum value to be 50, i.e., a GW event with SNR greater than 50 has the same color as the GW event with SNR of 50. In fact, many red dots have SNRs greater than 50, e.g., SNR at z = 0.146 for ET is 105.9, while for CE is 158.2. The reason of the design is that we can better highlight the comparison between ET and CE from the figure. If we set the maximal value of SNR higher (to about 100), almost all the data points are in blue, so it is of d L from CE are smaller than those from ET; (ii) SNR of CE is higher than that of ET at the same redshift. In Fig. 3, we show the simulated GW standard sirens from ET and CE. In the left panel, we show the standard siren data points without Gaussian randomness, where the central value of the luminosity distance is calculated by the fiducial cosmological model. In the right panel, in order to reflect the fluctuations in measured values resulting from actual observations, we show the standard siren data points with Gaussian randomization (the central values are populated according to a Gaussian distribution with mean being d L and standard deviation being σ dL ). In principle, the right panel is more representative of actual observational data, but the central values of d L have no effect on determining the absolute errors of cosmological parameters. Therefore, we use the data points in the left panel to constrain the cosmological models, because this is more helpful in investigating how the parameter degeneracies are broken to improve measurement precisions of cosmological parameters. We can clearly see that the measurement errors of d L from CE are smaller than those from ET, because CE has a better sensitivity than ET.

B. Other cosmological observations
In this work, we consider three current mainstream EM cosmological observations, including CMB, BAO, and SN. For the CMB data, we consider the Planck TT, TE, EE spectra at ≥ 30, the low-temperature Commander likelihood, and the low-SimAll EE likelihood from the Planck 2018 release [1]. For the BAO data, we consider the measurements from 6dFGS (z eff = 0.106) [199], SDSS-MGS (z eff = 0.15) [200], and BOSS DR12 (z eff = 0.38, 0.51, and 0.61) [201]. For the SN data, we use the latest Pantheon sample, which is comprised of 1048 data points from the Pantheon compilation [202].

C. Method of constraining cosmological parameters
To resolve the large-scale instability problem in the IDE cosmology [203], we apply the ePPF approach [182,183] for the IDE scenario so that the whole parameter space of IDE models can be explored without any divergence of the dark-energy perturbation. In this work, we employ the modified version of the available Markov-Chain Monte Carlo package CosmoMC [204], with the ePPF code [182,183] inserted, to constrain the neutrino mass and other cosmological parameters. In order to show the impacts of GW data difficult to visually compare SNRs of the two detectors from the figure. from CE and ET on constraining cosmological parameters, we use CMB+BAO+SN, CMB+BAO+SN+ET, and CMB+BAO+SN+CE to make our analysis. For convenience, we use CBS to standard for CMB+BAO+SN in the following. For the GW standard siren observation with N data points, the χ 2 function can be written as where z i , d i L , and σ i dL are the ith GW event's redshift, luminosity distance, and the measurement error of the luminosity distance, respectively. Ω represents the set of cosmological parameters.
When considering the combination of the current EM observations and the GW standard siren observations, the total χ 2 tot function is

III. RESULTS AND DISCUSSION
In this section, we report the constraint results of cosmological parameters in the ΛCDM+ m ν , IΛCDM+ m ν , and IwCDM+ m ν models. In these models, the three mass hierarchy cases of neutrinos, i.e., the NH, IH, and DH cases, have been considered. The constraint results of the NH case are shown as a representative in Figs. 4-6 and the constraint results are summarized in Tables I-III. Note that for the constraints on the total neutrino mass, the 2σ upper limits are given.
Note also that using the squared mass differences derived from the neutrino oscillation experiments, one can obtain the lower limits for the total neutrino mass, i.e., 0.05 eV for NH and 0.1 eV for IH; in the case of DH, the smallest value of the total neutrino mass is zero. For a parameter ξ, we use σ(ξ) and ε(ξ) to represent its absolute and relative errors, respectively, with ε(ξ) defined as ε(ξ) = σ(ξ)/ξ.
We first take a look at the results in the ΛCDM+ m ν model. In Fig. 4, we show the constraints on the ΛCDM+ m ν model in the m ν -H 0 and Ω m -H 0 planes from the CBS, CBS+ET, and CBS+CE data. We find that the addition of the GW data to the CBS data could lead to the reduction of the upper limits of m ν to some extent. The CBS+CE data give slightly smaller upper limits on m ν than those from the CBS+ET data. Concretely, when adding the ET data to the CBS data, the upper limits on m ν could be reduced by 2.7%-12.4% in the three hierarchy cases. While for CE, the upper limits on m ν could be reduced by 4.3%-14.0% in the three hierarchy cases. Here the results of ET are consistent with the previous results in Ref. [138].
Although using the GW data could only slightly improve the limits on the neutrino mass, they can significantly help improve the constraints on other cosmological parameters. We find that the constraints on Ω m and H 0 could be improved by 29.0%-32.8% and 30.4%-34.7%, respectively, when adding the ET data to the CBS data, and by 40.3%-43.8% and 43.5%-46.9%, respectively, for the case of CE.
In Fig. 5, we show the constraints on the IΛCDM+ m ν model in the m ν -β and Ω m -H 0 planes from the CBS, CBS+ET, and CBS+CE data. We can clearly see that, when considering the interaction be- tween vacuum energy and dark matter, the improvement of the limits on m ν by adding GW data is rather not evident. In the case of ET, the improvement of the limit on m ν is only 0.7%-1.8%, and in the case of CE, the improvement is 1.8%-4.1%. Therefore, we find that compared with the standard ΛCDM model, in its interaction version, the IΛCDM model, the improvement of the lim-its on m ν by GW data from ET and CE becomes weaker. This is because the IΛCDM model considers an extra cosmological parameter β compared with the ΛCDM model, which will degenerate with other cosmological parameters when the CBS data are used to constrain the IΛCDM model. Hence, compared with the ΛCDM model, the addition of the GW data to the CBS    data for its interaction version leads to weaker improvement.
We also find that the constraints on the coupling parameter β can be improved by using the GW data to a certain extent. In the IΛCDM+ m ν model, the constraints on β are improved by 19.2%-20.8% and 22.3%-26.2%, respectively, when the GW data of ET and CE are added on the basis of the CBS case.
In Fig. 6, we show the constraints on the IwCDM+ m ν model in the m ν -w and w-β planes from the CBS, CBS+ET, and CBS+CE data. We find that in this case the improvement of the limits on the neutrino mass is better than in the previous case. For ET, the improvement of the limit on m ν is 2.0%-5.4%, and for CE, the improvement is 5.3%-8.7%.
We find that in this case the constraints on the coupling parameter β and the EoS parameter of dark energy w can both be significantly improved by considering the addition of GW data. The constraints on β and w are improved by 2.4%-8.0% and 10.8%-13.2%, respectively, when considering the ET data, and by 7.1%-10.2% and 18.9%-21.1%, respectively, when considering the CE data.
In this work, we discuss the cosmological constraints on the IDE models in the cases of considering the GW standard siren observations from 3G ground-based GW detectors ET and CE. The results show that the limits on the neutrino mass can only be slightly improved with the help of the GW data, on the basis of the CBS constraint. Since the GW data can precisely constrain the Hubble constant H 0 , the addition of them in the cosmological fit can help break the cosmological parameter degeneracies formed by other cosmological observations. Therefore, the consideration of GW standard siren data can help significantly improve the constraints on the most cosmological parameters. However, the effect of massive neutrinos in the late universe and on the large scales can-not be distinctively distinguished from that of the cold dark matter, leading to the improvement of the limits on the neutrino mass by considering GW data is not obvious. Anyway, even though the impact on constraining the neutrino mass is not apparent, the GW standard sirens are rather useful in helping improve the constraints on the most cosmological parameters including the EoS of dark energy and the coupling between dark energy and dark matter.

IV. CONCLUSION
In the era of 3G ground-based GW detectors, a lot of GW standard siren data with the know redshifts could be obtained by the multi-messenger observation for BNS merger events. Obviously, these standard sirens would exert great impacts on the cosmological parameter estimation. Since the GW standard sirens can tightly constrain the Hubble constant, the consideration of them in a joint cosmological fit can lead to the cosmological parameter degeneracies formed by other cosmological observations being well broken. The GW standard sirens can thus be used to help significantly improve the constraints on cosmological parameters in the future.
It is of great interest to investigate whether the limits on the total neutrino mass can also be effectively improved by considering the GW standard siren data. In particular, the cosmological constraints on the neutrino mass are strongly model-dependent, and so the cases in different cosmological scenarios are needed to be detailedly discussed. In this work, we discuss the issue of weighing neutrinos in the IDE models by using the GW standard siren observations by ET and CE.
We consider the simplest IDE models, namely the IΛCDM and IwCDM models with Q = βHρ c . We simulate the GW standard siren data of the BNS mergers observed by ET and CE (in a way of multi-messenger detection). We investigate whether the GW standard sirens observed by ET and CE could help improve the constraint on the neutrino mass in the IDE models.
It is found that the GW standard siren observations from ET and CE can only slightly improve the constraint on the neutrino mass in the IDE models, compared to the current limit given by CMB+BAO+SN. This is mainly because the effect of massive neutrinos in the late universe and on the large scales cannot be distinctively distinguished from that of the CDM, leading to the improvement of the limits on the neutrino mass by considering GW data is not obvious. Although the limit on neutrino mass can only be slightly updated by considering the GW standard sirens, they are fairly useful in helping improve the constraints on the most cosmological parameters in-cluding the EoS of dark energy and the coupling between dark energy and dark matter.