Prospects for weighing neutrinos in interacting dark energy models using joint observations of gravitational waves and γ-ray bursts

Experiments involving solar and atmospheric neutrinos have demonstrated that neutrinos possess mass and exhibit substantial mixing among different species [1]. Despite this, directly measuring the absolute mass scale of neutrinos remains a formidable challenge in particle physics experiments. Neutrino mass has implications for the cosmic microwave background (CMB) and the universe's large-scale structure [2], which makes cosmological observations a crucial approach for determining their absolute mass.

Current principal cosmological observations fall into two categories. Firstly, there are observations related to the universe's expansion history, such as baryon acoustic oscillations (BAO) and type Ia supernovae (SN). Secondly, observations concerning the cosmic structure's growth, such as redshift space distortions, weak gravitational lensing, and galaxy clusters counts, are also pivotal. These observations help constrain the total neutrino mass, as extensively discussed in Refs. [3−49].

Recently, Ref. [44] explored the cosmological constraints on total neutrino mass within the framework of vacuum energy interacting with cold dark matter (IΛCDM), employing current observations. That study suggested that considering a direct interaction might alter the upper limit of the total neutrino mass, indicating that interactions between vacuum and cold dark matter could influence the cosmological measurement of neutrino mass. Nevertheless, current cosmological observations do not yet tightly constrain the total neutrino mass [50, 51]. However, gravitational-wave (GW) standard siren observations show promise as a powerful cosmological probe to aid in constraining the total neutrino mass.

The era of multi-messenger astronomy was inaugurated by the observation of the binary neutrino star (BNS) merger event GW170817 and its electromagnetic (EM) counterparts [52, 53]. Analyzing the GW waveform from such events directly yields the absolute luminosity distance, known as a standard siren. Additionally, identifying the EM counterpart of a GW source allows for the measurement of its redshift. Establishing the relationship between cosmic distances and redshifts paves the way for probing the universe’s expansion history. Thus, GW standard sirens are expected to significantly influence the estimation of cosmological parameters.

In the coming decades, third-generation (3G) ground-based GW detectors, such as Cosmic Explorer (CE) [54] and Einstein Telescope (ET) [55], are set to revolutionize the measurement of GW signals. These detectors are anticipated to form a network that will significantly enhance the capability to detect GW events. Ref. [56] presented a preliminary discussion on using future GW standard siren observations to weigh neutrinos within interacting dark energy (IDE) models. This analysis is based on an assumption that 1000 standard sirens could be detected over a 10-year period by either the ET or CE. While this estimate aligns with projections in other studies [29, 43, 44], variations in the redshift distribution could impact cosmological parameter estimations.

This work aims to explore the potential of 3G standard siren observations to measure neutrino mass in IDE models. We will analyze the detection strategy employing a network comprising the ET and two CEs — one in the US with a 40-km arm length and another in Australia with a 20-km arm length, collectively referred to as ET2CE. Additionally, we consider observations of γ-ray bursts (GRBs) from a detector similar to THESEUS, enhancing our capability to ascertain cosmological parameters.

In our analysis of IDE models, we focus exclusively on the IΛCDM models, characterized by a vacuum energy equation of state parameter $ w = -1 $. In these models, the energy conservation equations for vacuum energy and cold dark matter are described as follows:

where $ \rho_{\rm{de}} $ and $ \rho_{\rm{c}} $ represent the densities of dark energy and cold dark matter, respectively, H is the Hubble parameter, the dot denotes the derivative with respect to cosmic time t, and Q is the rate of energy transfer. In the field of IDE, one common assumption is that Q is proportional to the density of either dark matter or dark energy, i.e., $ Q = \beta H \rho_{\rm{c}} $ or $ Q = \beta H \rho_{\rm{de}} $ [57−66], where β is a dimensionless coupling parameter and H facilitates computational ease. Another assumption posits Q as $ Q = \beta H_0 \rho_{\rm{c}} $ or $ Q = \beta H_0 \rho_{\rm{de}} $ [38, 40, 67−71], with the Hubble constant $ H_0 $ used for dimensional consistency. In this work, we consider two forms of the energy transfer rate: $ Q_1 = \beta H \rho_{\rm{c}} $ and $ Q_2 = \beta H_0 \rho_{\rm{c}} $. According to Eqs. (1) and (2), a positive β implies that cold dark matter decays into vacuum energy, a negative β suggests that vacuum energy decays into cold dark matter, and $ \beta = 0 $ indicates no interaction between them.

This paper is organized as follows. In Section II, we describe the methodology in our analysis. In Section III, we report the constraint results and make some relevant discussions. In Section IV, we present a conclusion of this work.

The parameter space vector of the IΛCDM model is $ \{\omega_b,\; \omega_c,\; \tau,\; 100\theta_{\rm{MC}},\; \beta,\; \ln (10^{10}A_s),\; n_s\} $, where $ \omega_b $ and $ \omega_c $ are the present density of baryons and cold dark matter, respectively, τ is the Thomson scattering optical depth due to reionization, $ \theta_{\rm{MC}} $ (multiplied by 100) is the radio between the comoving sound horizon and angular diameter distance at the decoupling epoch, β is the coupling parameter in the IΛCDM model, $ A_s $ is the amplitude of primordial scalar perturbation power spectrum, and $ n_s $ is its power-law spectral index. When neutrino mass is considered in the IΛCDM model, one extra free parameter $ \sum m_\nu $ should be involved in the calculation.

Because we add neutrino mass into the IΛCDM model, the model considered in this paper is called the IΛCDM+$ \sum m_\nu $ model. For convenience, in this paper, we use IΛCDM1+$ \sum m_\nu $ and IΛCDM2+$ \sum m_\nu $ to denote the corresponding $ Q_1 = \beta H \rho_{\rm{c}} $ and $ Q_2 = \beta H_0 \rho_{\rm{c}} $ models that involve neutrino mass, respectively. Thus, there are eight independent parameters in total for these IΛCDM+$ \sum m_\nu $ (IΛCDM1+$ \sum m_\nu $ and IΛCDM2+$ \sum m_\nu $) models.

It should be mentioned that, when we consider the situation of vacuum energy interacting with cold dark matter, then the vacuum energy is not a pure background any more, and in this case, we must consider the vacuum energy perturbations. Note that, in the calculation of the dark energy perturbation evolution, a problem of perturbation instability appears; namely, the cosmological perturbations of dark energy will be divergent in some parts of the parameter space, which ruins the IDE cosmology at the perturbation level. Under such a circumstance, to avoid the perturbation instability problem in the IΛCDM model, in this work, we treat the vacuum energy perturbations based on the extended parameterized post-Friedmann (PPF) approach [72−74] (for the original version of the PPF method, see Refs. [75, 76]). For more information about the calculation of cosmological perturbations in the IDE models and the PPF approach, see Refs. [29, 65, 77, 78].

In this paper, we simulate the GW standard siren data from the 3G GW detectors and future GRB detector and use them to constrain neutrino mass in the IΛCDM models. We will investigate whether the GW standard sirens can improve the constraint on neutrino mass.

To show the constraining capability of the simulated GW standard siren data, we consider two data combinations for comparison: CMB+BAO+SN (abbreviated as CBS) and CBS+GW. For the CMB data, we use the CMB likelihood, including the TT, TE, and EE spectra at $ l\geq 30 $, low-l temperature commander likelihood, and low-l SimAll EE likelihood, from the Planck 2018 data release [51]. For the BAO data, we use the measurements from the 6dFGS ($ z_{\rm{eff}} = 0.106 $) [79], SDSS-MG ($ z_{\rm{eff}} = 0.15 $) [80], and BOSS-DR12 ($ z_{\rm{eff}} = $ 0.38, 0.51, and 0.61) [81]. For the SN data, we use the Pantheon sample, which is comprised of 1048 data points [82].

For the GW data, we use the simulation method described in Ref. [83]. Here, we provide only a brief overview. It is important to note that our approach includes a comprehensive calculation of the redshift distribution of GW-GRB events. This differs from the method used in previous studies, which assumed the detection of 1000 bright sirens over a 10-year observation period, as discussed in Refs. [56, 84−96].

Based on the star formation rate [97−99], the merger rate in the observer frame is

where ${\rm d}V(z)/{\rm d}z$ is the comoving volume element. $ {\cal{R}}_{\rm{m}} $ is the BNS merger rate in the source frame that is given by

where $ {\cal{R}}_{\rm{f}} $ is the cosmic star formation rate, for which we adopt the Madau-Dickinson model of Ref. [100], $ t(z) $ is the universe's age at the time of merger, $ t_{\rm{d}} $ is the time delay, and $ t_{\rm{min}} $ and $ t_{\rm{max}} $ are the minimum and maximum delay time for a massive binary to evolve merger [97], respectively. Here, $ t_{\rm{min}} = 20 $ Myr and $ t_{\rm{max}} $ is the Hubble time. The overall normalization is fixed by requiring that the value of ${\cal{R}}_{\rm{m}}(z = 0) = 920\; {\rm{Gpc}}^{-3}\; {\rm{yr}}^{-1}$ agrees with the local rate estimated from the O1 LIGO and O2 LIGO/Virgo observation run [101], which is also in accordance with the latest O3 observation run [102]. $ P(t_{\rm{d}}) $ is the distribution of the time delay $ t_{\rm{d}} $, and we follow and adopt the exponential form [97]

with an e-fold time of $ \tau = 100 $ Myr for $ t_{\rm{d}}>t_{\rm{min}} = 20 $ Myr. We simulate a catalog of BNS merger for 10-year observation. For each simulated GW source, the location $ (\theta,\phi) $, polarization angle ψ, coalescence phase $ \psi_{\rm{c}} $, and cosine of the inclination angle ι are drawn from uniform distributions. The mass of neutron star (NS) is assumed to be a Gaussian distribution, and the center value of the NS mass is $ 1.33\; M_{\odot} $, with a standard deviation of $ 0.09\; M_{\odot} $, where $ M_{\odot} $ is the solar mass.

Under the stationary phase approximation (SPA) [103], the Fourier transform of the frequency-domain waveform for a detector network (with N detectors) can be written as [104, 105]

where Φ is the $ N\times N $ diagonal matrix with $\Phi_{ij} = 2\pi f\delta_{ij}({{\boldsymbol{n}}\cdot {\boldsymbol{r}}}_k)$, n is the GW propagation direction, and $ {\boldsymbol{r}}_k $ is the location of the k-th detector. h(f) is

where $S_{n,k}(f)$ is the one-sided noise power spectral density of the k-th detector, and $ h_{k}(f) $ is the frequency domain GW waveform.

We consider the waveform in the inspiralling stage for the non-spinning BNS system in this study, and we adopt the restricted Post-Newtonian (PN) approximation and calculate the waveform to the 3.5 PN order [106, 107]. The Fourier transform of the GW waveform of the k-th detector can be expressed as

where the detailed forms of $ \Psi(f/2) $ and $ \varphi_{k,(2,0)} $ can be found in Refs. [105, 106]. $ {\cal{A}}_k $ is the Fourier amplitude, which can be written as

where $ d_{\rm{L}} $ is the luminosity distance of the GW source, $ F_{+,k} $ and $ F_{\times,k} $ are the antenna response functions of the k-th GW detector, $ {\cal{M}}_{\rm{chirp}} = (1+z)\eta^{3/5}M $ is the observed chirp mass, $ M = m_1+m_2 $ is the total mass of binary system with component masses $ m_1 $ and $ m_2 $, and $ \eta = m_1 m_2/M^2 $ is the symmetric mass ratio. In this study, we adopt the GW waveform in the frequency domain, and then the time t is replaced by $ t_{f} = t_{\rm{c}}-(5 / 256) {\cal{M}}_{\rm{chirp}}^{-5 / 3}(\pi f)^{-8 / 3} $ [105, 106], where $ t_{\rm{c}} $ is the coalescence time.

To obtain the detection of GW events, we need to calculate the signal-to-noise ratio (SNR) for each GW event. The combined SNR for the detection network of N detectors is given by

The inner product is defined as

where $ f_{\rm{lower}} $ is the lower cutoff frequency ($ f_{\rm{lower}} = 1 $ Hz for ET and $ f_{\rm{lower}} = 5 $ Hz for CE), $ f_{\rm{upper}} = 2/(6^{3/2}2\pi M_{\rm{obs}}) $ is the frequency at the last stable orbit with $M_{\rm{obs}} = (m_1+ m_2)(1+z)$, a and b are column matrices of the same dimension, and the star represents the conjugate transpose operator. Here, we adopt the SNR threshold to be 12 in the simulation, and we assume a running period of 10 years and duty cycle of 100% for the GW detector network.

For a GRB detected in coincidence with a GW signal, we require the peak flux to be higher than the flux limit of the satellite. Based on the fitting results of GRB170817A [108], we adopt the Gassian structured jet profile model

where $ L_{\rm{iso}}(\theta_{\rm{v}}) $ is the isotropically equivalent luminosity of short GRB observed at different viewing angle $ \theta_{\rm{v}} $, and we assume the jet's direction is aligned with the binary orbital angular momentum, namely, $ \iota = \theta_{\rm{v}} $. $ L_{\rm{on}} $ is the on-axis isotropic luminosity defined by $ L_{\rm{on}} = L_{\rm{iso}}(0) $, and $ \theta_{\rm{c}} = 4.7^{\circ} $ is the characteristic angle of the core. The detection probability of a short GRB is determined by $ \Phi(L){\rm{d}}L $. $ \Phi(L) $ is the intrinsic luminosity function, given by

where L is the isotropic rest frame peak luminosity in the $ 1-10000 $ keV energy range, $ L_{*} $ is a characteristic value separating the two regimes, and α and β are the characteristic slopes describing these regimes. In this paper, we use the values $ L_{*} = 2\times10^{52} $ erg/sec, $ \alpha = -1.95 $, and $ \beta = -3 $, which is in accordance with Ref. [109]. In addition, we also assume a standard low-end cutoff in luminosity of $ L_{\rm{min}} = 10^{49} $ erg/sec. Here, we term the on-axis isotropic luminosity $ L_{\rm{on}} $ as the peak luminosity L. According to the relation between luminosity and flux for GRB [110, 111], we can convert the flux limit $ P_{\rm{T}} $ to the luminosity $ L_{\rm{iso}} $ [83], where $P_{\rm{T}} = 0.2$ ph/(s·cm²) is the 50−300 keV band for THESEUS-like telescope [112]. For the THESEUS-like telescope, we make the assumption of an 80% duty cycle. From the GW catalogue which has SNR larger than the threshold 12, then we can calculate the probability of the GRB detection for every GW event according to the probability distribution $ \Phi(L){\rm{d}}L $.

For the total uncertainty of the luminosity distance $ d_{\rm{L}} $, we consider the instrument error $ \sigma_{d_{\rm{L}}}^{\rm{inst}} $, the weak-lensing error $ \sigma_{d_{\rm{L}}}^{\rm{lens}} $, and the peculiar velocity error $ \sigma_{d_{\rm{L}}}^{\rm{pv}} $, and thus the total error is given by

We first estimate the instrument error $ \sigma_{d_{\rm{L}}}^{\rm{inst}} $ by using the Fisher information matrix. The Fisher information matrix of a network including N independent GW detectors can be written as

where $ \theta_i $ denotes nine parameters ($ d_{\rm{L}} $, $ {\cal{M}}_{\rm{chirp}} $, η, θ, ϕ, ι, $ t_{\rm{c}} $, $ \psi_{\rm{c}} $, ψ) for a GW event. The instrumental error of the parameter $ \theta_i $ can be calculated by $ \Delta \theta_i = \sqrt{(F^{-1})_{ii}} $, where $ F_{ij} $ is the total Fisher information matrix for the network of N interferometers.

The error caused by weak lensing is adopted from Refs. [113, 114]

The error caused by the peculiar velocity of the GW source is given by [115]

where $ H(z) $ is the Hubble parameter, c is the speed of light in a vacuum, and $ \sqrt{\langle v^2\rangle} $ is the peculiar velocity of the GW source; we roughly use $\sqrt{\langle v^2\rangle} =$ 500 km/s.

In this work, we wish to study the constraint on neutrino mass in the IΛCDM+$ \sum m_\nu $ models from GW standard sirens. As the first attempt to explore the impacts of GW-GRB joint observation on the neutrino mass, we do not consider all the different cases of 3G GW observations, i.e., single ET, single CE, the CE-CE network (one in the US and another in Australia), and the ET-CE-CE network (one ET detector and two CE-like detectors). Instead, we only consider one case of ET-CE-CE network (abbreviated as GW hereafter) as a concrete example to analyze.

In the following, we use these observational data to place constraints on these IΛCDM+$ \sum m_\nu $ models. We will use two data combinations, i.e., CBS and CBS+GW, to constrain the neutrino mass $ \sum m_\nu $ and other cosmological parameters. We employ the modified version of the publicly available Markov-Chain Monte Corlo (MCMC) package CosmoMC [116] to perform the calculations, and we use the PPF package [72, 73] for the IΛCDM models to handle the perturbations of vacuum energy.

In this section, we report the fitting results of these IΛCDM+$ \sum m_\nu $ models and discuss the impacts of GW standard siren observation on constraining neutrino mass. We found that only approximately 640 standard sirens can be detected for the ET2CE+THESEUS network in 10 years, whose redshift distribution is given in Ref. [83]. The fitting results are presented in Tables 1 and 2 as well as Figs. 1 and 2. In Table 1, we quote $ \pm 1\sigma $ errors for the parameters, but for the parameters that cannot be well-constrained, e.g., neutrino mass $ \sum m_\nu $, we quote the 95.4% CL upper limits.

Parameter	ΛCDM+$ \sum m_\nu $ ($ Q = 0 $)			IΛCDM1+$ \sum m_\nu $ ($ Q_1 = \beta H\rho_{\rm{c}} $)			IΛCDM2+$ \sum m_\nu $ ($ Q_2 = \beta H_0\rho_{\rm{c}} $)
	CBS	CBS+GW		CBS	CBS+GW		CBS	CBS+GW
$ \Omega_m $	$ 0.3097^{+0.0062}_{-0.0067} $	$ 0.3098\pm0.0012 $		$ 0.3078\pm0.0081 $	$ 0.3076\pm0.0013 $		$ 0.3030^{+0.0160}_{-0.0180} $	$ 0.3024^{+0.0048}_{-0.0041} $
$ \sigma_8 $	$ 0.8124^{+0.0128}_{-0.0088} $	$ 0.8136^{+0.0117}_{-0.0086} $		$ 0.8140^{+0.0140}_{-0.0130} $	$ 0.8150^{+0.0130}_{-0.0110} $		$ 0.8190\pm0.0200 $	$ 0.8200\pm0.0140 $
$ H_0 $	$ 67.75^{+0.52}_{-0.48} $	$ 67.75\pm0.05 $		$ 67.92\pm0.65 $	$ 67.92\pm0.05 $		$ 68.06^{+0.81}_{-0.82} $	$ 68.06\pm0.06 $
β	...	...		$ 0.0005\pm0.0013 $	$ 0.0005^{+0.0009}_{-0.0011} $		$ 0.0230^{+0.0470}_{-0.0480} $	$ 0.0230^{+0.0260}_{-0.0290} $
$ \sum m_\nu $	$<0.123 $	$<0.103 $		$<0.150 $	$<0.143 $		$<0.156 $	$<0.140 $

Table 1.  Fitting results for the IΛCDM+$ \sum m_\nu $ models using the CBS and CBS+GW data. We quote $ \pm 1\sigma $ errors, but for the parameters that cannot be well-constrained, we quote the 95.4% CL upper limits. Here, $ H_0 $ is in units of $ {\rm{km}}\ { {\rm{s}}^{-1}}\ { {\rm{Mpc}}^{-1}} $, and $ \sum m_\nu $ is in units of eV.

Parameter	ΛCDM+$ \sum m_\nu $ ($ Q = 0 $)			IΛCDM1+$ \sum m_\nu $ ($ Q_1 = \beta H\rho_{\rm{c}} $)			IΛCDM2+$ \sum m_\nu $ ($ Q_2 = \beta H_0\rho_{\rm{c}} $)
	CBS	CBS+GW		CBS	CBS+GW		CBS	CBS+GW
$ \sigma(\Omega_m) $	$ 0.0065 $	$ 0.0012 $		$ 0.0081 $	$ 0.0013 $		$ 0.0170 $	$ 0.0045 $
$ \sigma(\sigma_8) $	$ 0.0108 $	$ 0.0102 $		$ 0.0135 $	$ 0.0120 $		$ 0.0200 $	$ 0.0140 $
$ \sigma(H_0) $	$ 0.50 $	$ 0.05 $		$ 0.65 $	$ 0.05 $		$ 0.82 $	$ 0.06 $
$ \sigma(\beta) $	$ ... $	$ ... $		$ 0.0013 $	$ 0.0010 $		$ 0.0475 $	$ 0.0275 $
$ \varepsilon(\Omega_m) $	$ 2.10$%	$ 0.39$%		$ 2.63$%	$ 0.42$%		$ 5.61$%	$ 1.49$%
$ \varepsilon(\sigma_8) $	$ 1.33$%	$ 1.25$%		$ 1.66$%	$ 1.47$%		$ 2.44$%	$ 1.71$%
$ \varepsilon(H_0) $	$ 0.74$%	$ 0.07$%		$ 0.96$%	$ 0.07$%		$ 1.20$%	$ 0.09$%

Table 2.  Absolute and relative errors of cosmological parameters in the ΛCDM+$ \sum m_\nu $ and IΛCDM+$ \sum m_\nu $ models using the CBS and CBS+GW data combinations. Here, $ H_0 $ is in units of $ {\rm{km}}\ { {\rm{s}}^{-1}}\ { {\rm{Mpc}}^{-1}} $. For a parameter ξ, $ \sigma(\xi) $ and $ \varepsilon(\xi) = \sigma(\xi)/\xi $ represent its absolute and relative errors, respectively.

Figure 1.  (color online) Two-dimensional marginalized posterior contours (68.3% and 95.4% CL) in the $ \Omega_{\rm{m}} $–$ \sum m_\nu $ plane for the IΛCDM+$ \sum m_\nu $ models using the CBS (left) and CBS+GW (right) data combinations.

Figure 2.  (color online) Two-dimensional marginalized posterior contours (68.3% and 95.4% CL) in the $ \sum m_\nu $–β and $ \sum m_\nu $–$ H_0 $ planes for the IΛCDM1+$ \sum m_\nu $($ Q_1 = \beta H\rho_{\rm{c}} $) (left) and IΛCDM2+$ \sum m_\nu $($ Q_2 = \beta H_0\rho_{\rm{c}} $) (right) models using the CBS and CBS+GW data combinations.

First, we discuss the effect of simulated GW data of ET2CE on constraining the total neutrino mass. In the case of CBS constraints, from Table 1, we have $ \sum m_\nu<0.123 $ eV for ΛCDM+$ \sum m_\nu $, $ \sum m_\nu<0.150 $ eV for IΛCDM1+$ \sum m_\nu $, and $ \sum m_\nu<0.156 $ eV for IΛCDM2+$ \sum m_\nu $. After adding the GW data of ET2CE, namely, using the data combination CBS+GW, the constraint results become $ \sum m_\nu<0.103 $ eV for ΛCDM+$ \sum m_\nu $, $ \sum m_\nu<0.143 $ eV for IΛCDM1+$ \sum m_\nu $, and $ \sum m_\nu<0.140 $ eV for IΛCDM2+$ \sum m_\nu $. We found that, compared with the ΛCDM+$ \sum m_\nu $ model, the constraint on neutrino mass $ \sum m_\nu $ becomes slighter looser in these IΛCDM+$ \sum m_\nu $ models. We show the two-dimensional marginalized posterior contours (68.3% and 95.4% CL) in the $ \Omega_{\rm{m}} $–$ \sum m_\nu $ plane for the ΛCDM+$ \sum m_\nu $and IΛCDM+$ \sum m_\nu $ models in Fig. 1 (left panel for CBS and right panel for CBS+GW). From this figure, we can clearly see that once the interaction is considered in the model, the upper limit of $ \sum m_\nu $ becomes larger.

Moreover, we also found that the GW data help reduce the upper limits of neutrino mass $ \sum m_\nu $ by 16.26%, 4.67%, and 10.26% for the ΛCDM+$ \sum m_\nu $, IΛCDM1+$ \sum m_\nu $, and IΛCDM2+$ \sum m_\nu $ models, respectively. When the GW data are considered, the constraints on $ \sum m_\nu $ become slightly tighter in the ΛCDM+$ \sum m_\nu $ and IΛCDM+$ \sum m_\nu $ models (also see Fig. 2). Thus, in the IΛCDM+$ \sum m_\nu $ models, compared to the current observations, the GW standard siren observations from ET2CE can slightly improve the constraint on neutrino mass, which is in accordance with the conclusion on the IDE models in the previous study [56].

Next, we discuss how the GW data help improve the constraints on other parameters, i.e., $ \Omega_m $, $ \sigma_8 $, and $ H_0 $. In Table 2, we show the absolute and relative errors of $ \Omega_m $, $ \sigma_8 $, and $ H_0 $ from the CBS and CBS+GW data combinations. Comparing the results from two data combinations, we found that the constraints on $ \Omega_m $, $ \sigma_8 $, and $ H_0 $ are improved by 73.53%–83.95%, 5.56%–30%, and 90%– 92.68%, respectively, when considering the ET2CE data. Obviously, the GW data can indeed effectively improve the constraints on the parameters $ \Omega_m $ and $ H_0 $ (also see Figs. 1 and 2).

Finally, we present the constraint results of the coupling constant β. By using CBS data, we have $ \beta = 0.0005\pm0.0013 $ for the IΛCDM1+$ \sum m_\nu $ model and $ \beta = 0.0230^{+0.0470}_{-0.0480} $ for IΛCDM2+$ \sum m_\nu $ model. It is shown that a positive value of β is slightly favored and $ \beta>0 $ is only at $ 0.38\sigma $ and $ 0.48\sigma $ levels, respectively. By using CBS+GW data, we have $ \beta = 0.0005^{+0.0009}_{-0.0011} $ for the IΛCDM1+ $ \sum m_\nu $ model and $ \beta = 0.0230^{+0.0260}_{-0.0290} $ for the IΛCDM2+$ \sum m_\nu $ model. Thus, we found that the error is slightly shrunk and now $ \beta>0 $ is favored at $ 0.45\sigma $ and $ 0.79\sigma $ levels, respectively. For these two IΛCDM+$ \sum m_\nu $ models, one can clearly see that, in the IΛCDM+$ \sum m_\nu $ cases (with $ Q_1 = \beta H\rho_{\rm{c}} $ and $ Q_2 = \beta H_0\rho_{\rm{c}} $), $ \beta = 0 $ is favored within the $ 1\sigma $ significance range, no matter whether the GW data are taken into account. Thus, in the two IΛCDM+$ \sum m_\nu $ models, there is no evidence of a nonzero interaction. Moreover, comparing the values of β from the two data combinations, we also found that the accuracy of β is increased by 23.08% ($ Q_1 = \beta H \rho_{\rm{c}} $) and 42.11% ($ Q_2 = \beta H_0\rho_{\rm{c}} $) when adding the GW data of ET2CE in the cosmological fit (also see Fig. 2). This indicates that the GW data of the ET2CE can also improve the constraint accuracies of coupling parameter β.

In this paper, we investigated how GW standard sirens from 3G ground-based GW detectors can constrain the total neutrino mass ($ \sum m_\nu $) in interacting vacuum energy models with energy transfer forms $ Q_1 = \beta H\rho_{\rm{c}} $ and $ Q_2 = \beta H_0\rho_{\rm{c}} $. We used an extended version of the PPF approach to compute the vacuum energy perturbations within these interacting scenarios. We focused on the synergy between 3G GW detectors and a GRB detector for multi-messenger observations, specifically employing the ET2CE network as a representative for GW discussions. To evaluate the impact of GW data on the constraints of $ \sum m_\nu $, we also considered existing CMB+BAO+SN data for comparison and combination.

Our analysis revealed that GW data can reduce the upper limits of $ \sum m_\nu $ by 16.26%, 4.67%, and 10.26% for the ΛCDM+$ \sum m_\nu $ , IΛCDM1+$ \sum m_\nu $, and IΛCDM2+$ \sum m_\nu $ models, respectively. Thus, GW standard siren data from ET2CE offer a modest improvement in constraining $ \sum m_\nu $ compared to CBS alone. For the derived parameters $ \Omega_m $ and $ H_0 $, incorporating GW data into the cosmological fit substantially enhances the precision: the accuracy of $ \Omega_m $ improved by 73.53%–83.95% and that of $ H_0 $ by 90%–92.68%. These significant enhancements underscore the value of including GW data from ET2CE. Moreover, based on the combined constraints from CBS and CBS+GW data, we found that GW data from ET2CE also significantly refine the accuracy of the coupling strength β.

We thank Shang-Jie Jin for helpful discussions.