Impacts of gravitational-wave standard siren observation of the Einstein Telescope on weighing neutrinos in cosmology

We investigate the impacts of the gravitational-wave (GW) standard siren observation of the Einstein Telescope (ET) on constraining the total neutrino mass. We simulate 1000 GW events that would be observed by the ET in its 10-year observation by taking the standard $\Lambda$CDM cosmology as a fiducial model. We combine the simulated GW data with other cosmological observations including cosmic microwave background (CMB), baryon acoustic oscillations (BAO), and type Ia supernovae (SN). We consider three mass hierarchy cases for the neutrino mass, i.e., normal hierarchy (NH), inverted hierarchy (IH), and degenerate hierarchy (DH). Using Planck+BAO+SN, we obtain $\sum m_\nu<0.175$ eV for the NH case, $\sum m_\nu<0.200$ eV for the IH case, and $\sum m_\nu<0.136$ eV for the DH case. After considering the GW data, i.e., using Planck+BAO+SN+GW, the constraint results become $\sum m_\nu<0.151$ eV for the NH case, $\sum m_\nu<0.185$ eV for the IH case, and $\sum m_\nu<0.122$ eV for the DH case. We find that the GW data can help reduce the upper limits of $\sum m_\nu$ by 13.7%, 7.5%, and 10.3% for the NH, IH, and DH cases, respectively. In addition, we find that the GW data can also help break the degeneracies between $\sum m_{\nu}$ and other parameters. We show that the GW data of the ET could greatly improve the constraint accuracies of cosmological parameters.


I. INTRODUCTION
On 17 August 2017, the signal of a gravitational wave (GW) produced by the merger of a binary neutronstar system (BNS) was detected for the first time [1], and the electromagnetic (EM) signals generated by the same transient source were also observed subsequently, which indicates that the age of gravitational-wave multimessenger astronomy is coming. The measurement of GWs from the merger of a binary compact-object system involves the information of absolute luminosity distance of the transient source [2], and thus if we can simultaneously accurately measure the GW and EM signals from the same merger event of a BNS or a binary system consisting of a neutron star (NS) and a black hole (BH), then we are able to establish a true luminosity distanceredshift (d L -z) relation. Therefore, the GW observations can serve as a cosmic "standard siren", which can be developed to be a new cosmological probe if we can accurately observe a large number of merger events of this class.
The current mainstream cosmological probes include the measurements of cosmic microwave background (CMB) anisotropies (temperature and polarization), baryon acoustic oscillations (BAO), type Ia supernovae (SNIa), and the Hubble constant, etc. In addition, there are also some observations for the growth history of largescale structure (LSS), such as the shear measurement of the weak gravitational lensing, the galaxy clusters number counts in light of Sunyaev-Zeldovich (SZ) effect, and the CMB lensing measurement, etc. When using these observational data to make cosmological parameter esti- * Corresponding author † Electronic address: zhangxin@mail.neu.edu.cn mation, some problems occur including mainly: (i) there are degeneracies between some parameters, and in some of which the correlations are rather strong, and (ii) there are apparent tensions between some observations. The GW standard siren observations have some peculiar advantages in breaking the parameter degeneracies, owing to the fact that the GW observations can directly measure the true luminosity distances, but the SNIa observations actually can only measure the ratios of luminosity distances at different redshifts, but not the true luminosity distances. In addition, compared to the standard candle provided by the SNIa observations [3][4][5], which needs to cross-calibrate the distance indicators on different scales, the GW observations allow us to directly measure the luminosity distances up to higher redshifts [6]. To obtain the information of redshifts, one needs to detect the EM counterparts of the GW sources. In fact, there are some forthcoming large facility observation programs, such as Large Synoptic Survey Telescope (LSST) [7], Square Kilometer Array (SKA) [8], and Extremely Large Telescope (ELT) [9], which can detect the EM counterparts by optically identifying the host galaxies. In this way, in the near future, the d L -z relation would be obtained by the GW standard siren measurements, and we could use this powerful tool to explore the expansion history of the universe. Recently, the related issues have been discussed by some authors [10][11][12][13][14][15][16][17][18][19] (see also Ref. [20] for a recent review). For example, in Ref. [10], Cai and Yang estimated the ability of using GW data to constrain cosmological parameters. They considered to use the GW detector under planning, the Einstein Telescope (ET), to simulate the GW data, which is a third-generation ground-based detector of GWs [21]. ET is ten times more sensitive than the current advanced ground-based detectors and it covers the frequency range of 1 − 10 4 Hz. According to their results [10], the errors of cosmological parameters can be constrained to be ∆h ∼ 5 × 10 −3 and ∆Ω m ∼ 0.02 when using 1000 GW events, whose sensitivity is comparable to that of the Planck data [22]. From their study, we can see that the GW data can indeed be used to improve the constraint accuracies of parameters. In Ref. [23], we can see how the parameter degeneracies are broken by the GW observations in an efficient way. In this work, we investigate the issue of measuring the neutrino mass in light of cosmological observations and we will discuss what role the GW observations of the ET will play as a new cosmological probe in this study.
Since the phenomenon of neutrino oscillation was discovered, which proved that the neutrino masses are not zero [24], the determination of neutrino masses has been an important issue in the field of particle physics. Due to the fact that neutrino oscillation experiments are only sensitive to the squared mass differences between the neutrino mass eigenstates, it is a great challenge to determine the absolute masses of neutrinos by particle physics experiments. Although the solar and atmospheric neutrino oscillation experiments can give two squared mass differences between the mass eigenstates: ∆m 2 21 7.5 × 10 −5 eV 2 and |∆m 2 32 | 2.5 × 10 −3 eV 2 , we cannot determine whether the third neutrino is heaviest or lightest. Therefore, these measurements can only give two possible mass orders, i.e., the normal hierarchy (NH) with m 1 < m 2 m 3 and the inverted hierarchy (IH) with m 3 m 1 < m 2 . If the total mass of neutrinos ( m ν ) can be measured, then the absolute masses of neutrinos could be solved by combining the total mass and the two squared mass differences.
Massive neutrinos play an significant role in the evolution of the universe, and thus they leave distinct signatures on CMB and LSS at different epochs of the evolution of the universe. These signatures can actually be extracted from the cosmological observations, from which the total mass of neutrinos can be effectively constrained. In recent years, the combinations of various cosmological observations have been providing more and more tight constraint limits for the total mass of neutrinos. For the latest progresses on this aspect, see e.g. Refs. [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43].
In fact, in a recent paper [44], the constraints on the total neutrino mass have been discussed by considering the inclusion of the actual observation of GW and EM emission produced by the merger of BNS, GW170817. In Ref. [44], the authors considered a 12-parameter extended cosmological model that contains the dark-energy equation-of-state parameter w 0 and w a , the spatial curvature Ω k , the total neutrino mass m ν , the effective number of relativistic species N eff , and the running of the scalar spectral index dn s /d ln k, besides the 6 base parameters, and they made a comparison for the constraint results of Planck and Planck+GW170817. For such a 12parameter model, using only Planck data gives m ν < 1.11 eV, and the combination of Planck+GW170817 gives m ν < 0.77 eV, showing that the inclusion of GW170817 leads to an about 30% improvement for the upper limit of m ν , compared to the result obtained from the Planck data alone. However, it is clearly known that the 12-parameter model actually cannot be well constrained by the current observations. In the present work, we wish to scrutinize the standard cosmology, i.e., a 7-parameter Λ cold dark matter (ΛCDM) model that contains the 6 base parameters plus m ν . We aim to see how the future GW observations help improve the constraints on the total neutrino mass in a 7-parameter ΛCDM cosmology.
In this study, our main focus is on the ability of the GW observations of ET to constrain the neutrino mass. We follow Ref. [10] to generate the catalogue of GW events, from which we can get the corresponding d L -z relation. According to the influences of instruments and weak lensing, we can estimate the uncertainty on the measurement of d L . We simulate 1,000 GW events that could be observed by the ET in its 10-year observation. Then, we combine these simulated GW data with other current observations to constrain the total mass of neutrinos, m ν , by using the Markov-chain Monte Carlo (MCMC) [45] approach.
The paper is organized as follows. In Sec. II, we describe the method to generate simulated GW events and get the luminosity distances with the simulated measurement errors. In addition, the fiducial model and the data processing method are also briefly introduced. In Sec. III, we report the results of this work. In Sec. IV, we make a conclusion for this work.

II. METHODS OF SIMULATING DATA AND CONSTRAINING PARAMETERS
A. Method of simulating data The first step for generating GW data is to simulate the redshift distribution of the sources. Following Refs. [10,16], the distribution takes the form where d C (z) is the comoving distance at redshift z and R(z) denotes the time evolution of the burst rate and takes the form [10, 46, 47] 3 through the fiducial model, i.e., the ΛCDM model. If we consider a flat Friedmann-Robertson-Walker universe, the Hubble parameter H(z) for the ΛCDM cosmology can be written as where Ω m and Ω r represent the fractional energy densities of matter and radiation, respectively. In the late universe, Ω r can be neglected, and where Ω b , Ω c , and Ω ν represent the fractional energy densities of baryons, cold dark matter, and neutrinos, respectively. Note that Ω ν can be expressed as where h is the reduced Hubble constant (the Hubble constant H 0 = 100h km s −1 Mpc −1 ). Since the radiation component can be ignored in the late universe, we set Ω r = 0 here. The luminosity distance d L can be calculated by where E(z) ≡ H(z)/H 0 . According to the redshift distribution of the GW sources, we can then use Eq. (5) to generate a catalogue of the GW sources. That is to say, the d L -z relation can be obtained for every GW event in a ΛCDM cosmology. The next step is to get the error of d L of the GW source, which is denoted by σ d L in this paper. We first need to generate the waveform of GWs. Because the GW amplitude depends on the luminosity distance d L , we can get the information of d L and σ d L from the amplitude of waveform.
Following Ref. [10], the strain in GW interferometers can be written as where ψ is the polarization angle and (θ,φ) are angles describing the location of the source relative to the detector. Here the antenna pattern functions of the ET (i.e. F + and F × ) are [16] F (1) For the other two interferometers, their antenna pattern functions can be derived from above equations because the interferometers have 60 • with each other.
Following Refs. [16,17], we can compute the Fourier transform H(f ) of the time domain waveform h(t), where A is the Fourier amplitude that is given by where M c = M η 3/5 is the "chirp mass", M = m 1 + m 2 is the total mass of coalescing binary with component masses m 1 and m 2 , and η = m 1 m 2 /M 2 is the symmetric mass ratio. Note that all the masses here refer to the observed mass rather than the intrinsic mass. The observed mass is larger than the intrinsic mass by a factor of (1 + z), i.e., M obs = (1 + z)M int . ι is the angle of inclination of the binary's orbital angular momentum with the line of sight. Due to the fact that the short gamma ray bursts (SGRBs) are expected to be strongly beamed, the coincidence observations of SGRBs imply that the binaries should be orientated nearly face on (i.e., ι 0) and the maximal inclination is about ι = 20 • . In fact, averaging the Fisher matrix over the inclination ι and the polarization ψ with the constraint ι < 90 • is roughly the same as taking ι = 0 in the simulation [17]. Thus, when we simulate the GW source we can take ι = 0. But, when we estimate the practical uncertainty of the measurement of d L , the impacts of the uncertainty of inclination should be taken into account. Actually, the consideration of the maximal effect of the inclination (between ι = 0 and ι = 90 • ) on the signal-to-noise (SNR) leads to a factor of 2 [see Eq. (13)]. The definitions of other parameters and the values of the parameters can be found in Ref. [10]. After knowing the waveform of GWs, we can then calculate the SNR. A GW detection is confirmed if it produces a combined SNR of at least 8 in ET [21,48] (see also Refs. [10,11,14,16]). The combined SNR for the network of three independent interferometers is where ρ (i) = H (i) , H (i) , with the inner product being defined as where a "∼" above a function denotes the Fourier transform of the function and S h (f ) is the one-side noise power spectral density. In this work, we take S h (f ) of the ET to be the same as in Ref. [16].
Using the Fisher information matrix, we can estimate the instrumental error on the measurement of d L , which can be written as Because we only focus on the parameter d L in the waveform, we find that σ inst d L d L /ρ due to H ∝ d −1 L . Considering the effect from the inclination angle ι, we add a factor 2 in front of the error, so the error is written as Following Ref. [10], we set the additional error σ lens d L = 0.05zd L , which represents the error from weak lensing. Thus, the total error on the measurement of d L can be expressed as Using the method described above, we can generate the catalogue of the GW events with their z, d L , and σ d L . According to the results in Ref. [10], we know that it requires more than 1000 GW events to match the Planck sensitivity, and so we simulate 1000 GW events that are expected to be detected by the ET in its 10-year observation.

B. Method of constraining parameters
In order to constrain cosmological parameters, we use the MCMC method to infer the posterior probability distributions of parameters. To measure the total neutrino mass in light of cosmological observations, we set m ν to be a free parameter in the cosmological fit.
In this work, we add the simulated GW data in the combined cosmological data. For the GW standard siren measurement with N simulated data points, we can write its χ 2 as wherez i ,d i L , andσ i d L are the ith redshift, luminosity distance, and error of luminosity distance of the simulated GW data, and Ω represents the set of cosmological parameters.
To show the constraining capability of the simulated GW data, we consider two data combinations for comparison in this work: (i) Planck+BAO+SN and (ii) Planck+BAO+SN+GW. For the CMB data, we use the Planck 2015 temperature and polarization data. For the BAO data, we use the measurements of the six-degreefield galaxy (6dFGS) at z eff = 0.106 [49], the SDSS main galaxy sample (MGS) at z eff = 0.15 [50], the baryon oscillation spectroscopic survey (BOSS) LOWZ at z eff = 0.32 [51], and the BOSS CMASS at z eff = 0.57 [51]. For the SN data, we use the "joint light-curve analysis" (JLA) sample [52]. For the simulated GW data, we consider 1,000 GW events that could be observed by the ET in its 10-year observation.
For the neutrino mass measurement in this work, we consider three mass hierarchy cases, i.e., the normal hierarchy (NH), the inverted hierarchy (IH), and the degenerate hierarchy (DH). For the details of this aspect, see Refs. [29,30].

III. RESULTS AND DISCUSSION
Our constraint results for the neutrino mass and other cosmological parameters are shown in Table I and Figs. 1-4. Note that in this work we have considered three mass hierarchy cases for massive neutrinos and we have used two data combinations to make the analysis.
For convenience, we also use "Data1" to represent the Planck+BAO+SN data combination, and use "Data2" to represent the Planck+BAO+SN+GW data combination; e.g., see Table I. In Table I, we show the best-fit results with the 68% CL uncertainties for the cosmological parameters, but owing to the fact that the neutrino mass cannot be well constrained, we only give the 95% CL upper limits for the neutrino mass m ν . In addition, the derived parameters Ω m and H 0 are also listed in this table.
In Fig. 1, we show the one-dimensional posterior distributions of m ν using the two data combinations, for the three mass hierarchy cases. In Figs. 2-4, we show the two-dimensional posterior distribution contours (68% and 95% CL) in the Ω m -H 0 , H 0 -m ν , and Ω mm ν planes, for the three mass hierarchy cases, also using the two data combinations. In these figures, the blue lines and contours represent the results from the Planck+BAO+SN data, and the red lines and contours represent the results from the Planck+BAO+SN+GW data.
First, we discuss the effect of the simulated GW data of the ET on constraining the total neutrino mass. From the one-dimensional posterior distributions of m ν in Fig. 1, we find that, when the GW data are considered, the constraints on m ν in all the three mass hierarchy cases become tighter, i.e., smaller values of the upper limits of m ν are obtained. The detailed constraint results have been given in Table I. Using the data combination Planck+BAO+SN, we obtain: m ν < 0.175 eV for the NH case, m ν < 0.200 eV for the IH case, and m ν < 0.136 eV for the DH case. After adding the GW data of the ET, namely when using the data combination Planck+BAO+SN+GW, the constraint results become: m ν < 0.151 eV for the NH case, m ν < 0.185 eV for   the IH case, and m ν < 0.122 eV for the DH case. We find that the GW data help reduce the upper limits of m ν by 13.7%, 7.5%, and 10.3% for the NH, IH, and DH cases, respectively. Obviously, the GW data can indeed effectively improve the constraints on the neutrino mass.
Next, we discuss how the GW data help improve the constraint accuracies for other cosmological parameters. The constraint results of the derived parameters Ω m and H 0 are also listed in Table I and H 0 = 67.92 ± 0.52 km s −1 Mpc −1 for the DH case. After considering the GW data of the ET, i.e., using the data combination of Planck+BAO+SN+GW, we obtain: Ω m = 0.3117 ± 0.0024 and H 0 = 67.52 ± 0.16 km s −1 Mpc −1 for the NH case, Ω m = 0.3158 ± 0.0026 and H 0 = 67.18 ± 0.17 km s −1 Mpc −1 for the IH case, and Ω m = 0.3090 ± 0.0023 and H 0 = 67.75 ± 0.14 km s −1 Mpc −1 for the DH case. Comparing the results from the two data combinations, we find that the accuracy of Ω m is increased by about 60% and the accuracy of H 0 is increased by about 68% when the GW data of the ET are considered in the cosmological fit. This indicates that the GW data of the ET can significantly improve the constraint accuracies of cosmological parameters.
We also display the two-dimensional posterior distribution contours in Figs. 2-4. From the blue contours (Planck+BAO+SN) in the H 0 -m ν and Ω m -m ν planes, we can see that m ν is anti-correlated with H 0 and is positively correlated with Ω m . From the red contours (Planck+BAO+SN+GW) in the H 0 -m ν and Ω m -m ν planes, we find that, when the GW data are considered, the parameter space is greatly shrunk in each plane and the constraints on Ω m and H 0 become much tighter. Moreover, after adding the GW data, there is no obvious correlation between m ν and other cosmological parameter. This indicates that the degeneracies be-tween parameters including the neutrino mass existing in the cosmological EM observations can be effectively broken by the GW observations. This is because when we use the Planck data to do the cosmological fit, the parameter combinations must be constrained to a constant θ * (the observed angular size of acoustic scale θ * = r s /D A ), which will cause the degeneracies between parameters. Also, neither BAO or SN can truly measure the cosmological distances (d A or d L ). But the GW data contain the absolute distance information at low redshifts relative to the CMB data, so they can help to break the degeneracies between m ν and other parameters. Hence, the upper limits on the total neutrino mass are also reduced.
To briefly summarize, our results show that the GW data can indeed improve the constraints on the total neutrino mass. When the GW data of the ET are considered, tighter bounds on the total neutrino mass could be obtained. Also, after considering the GW data, the constraints on the derived parameters Ω m and H 0 become much tighter. In addition, the GW data can help break the degeneracies between m ν and other parameters.

IV. CONCLUSION
In this paper, we investigated the constraint capability of the GW observation of the ET on the total neutrino mass m ν . We constrained the total neutrino mass in the ΛCDM cosmology by using the simulated GW data of ET in combination with other cosmological data including CMB, BAO and SN. For the three-generation neutrinos, we considered the cases of normal hierarchy, inverted hierarchy, and degenerate hierarchy. For the GW data, we considered the ΛCDM model as our fiducial model to simulate 1000 GW events that could be detected by the ET in its 10-year observation. In order to show the effect of the GW data, we used two data combinations including Planck+BAO+SN and Planck+BAO+SN+GW to constrain cosmological parameters.
Through comparing the constraint results from the two data combinations, we find that the GW data can indeed effectively improve the constraints on the total neutrino mass. Using Planck+BAO+SN, we obtain m ν < 0.175 eV for the NH case, m ν < 0.200 eV for the IH case, and m ν < 0.136 eV for the DH case. After considering the GW data, i.e., using Planck+BAO+SN+GW, the constraint results become m ν < 0.151 eV for the NH case, m ν < 0.185 eV for the IH case, and m ν < 0.122 eV for the DH case. It is found that the GW data can help reduce the upper limits of m ν by 13.7%, 7.5%, and 10.3% for the NH, IH, and DH cases, respectively. For the derived parameters Ω m and H 0 , the GW data can significantly improve the constraint accuracies of them. The accuracy of Ω m is increased by about 60% and the accuracy of H 0 is increased by about 68%, when considering the GW data in the cosmological fit. In addition, the GW data can help to break the degeneracies between m ν and other parameters.