Prospects for searching for sterile neutrinos in dynamical dark energy cosmologies using joint observations of gravitational waves and γ-ray bursts

On 17 August 2017, the signal of a gravitational wave (GW) produced by the merger of a binary neutron star system (BNS) was detected for the first time [1], and the electromagnetic (EM) signals generated by the same transient source were also observed subsequently, marking the beginning of the era of GW multi-messenger astronomy. GWs can act as "standard sirens" [2, 3] because the GW waveform encodes the absolute luminosity distance to the source. If the redshift of the EM counterpart of the GW source can be determined, a distance-redshift relation can be established, which can be used to probe the expansion history of the universe (see, e.g., Refs. [4−42] for related discussions).

In the future, third-generation (3G) ground-based GW detectors, such as the Einstein Telescope (ET) [43, 44] and the Cosmic Explorer (CE) [45, 46], are expected to detect a large number of standard sirens from BNS merger events. These observations can be used to constrain cosmological parameters with high precision, making GW standard sirens a powerful new cosmological probe. In particular, GW standard sirens can help break the parameter degeneracies generated by the current EM cosmological observations, thereby improving constraints on both active and sterile neutrino parameters; see, e.g., Refs. [12, 24, 36, 39].

The phenomenon of neutrino oscillation observed in solar and atmospheric neutrino experiments has shown that neutrinos are massive and that there is significant mixing among the three active neutrino flavors (see Ref. [47] for a review). However, anomalies observed in short-baseline neutrino oscillation experiments have hinted at the existence of additional light neutrino species [48−55], which, if they exist, would be sterile neutrinos. Explaining these anomalous results appears to require the existence of at least one sterile neutrino with mass at the eV scale. However, this result is not favored by several other neutrino experiments, such as the Daya Bay [56, 57] and IceCube [58, 59]. Consequently, further detailed investigations are necessary to explore the potential existence of sterile neutrinos. Meanwhile, laboratory experiments (such as JUNO, MINOS, and IceCube) provide relatively model-independent and direct tests but are limited by restricted sensitivity ranges, energy resolution, and statistical uncertainties [60−64]. By contrast, cosmological observations can yield much tighter constraints on sterile neutrino parameters through global fits to the evolution of the universe. Therefore, cosmological observations serve as an independent and complementary avenue for exploring the properties of sterile neutrinos. Extensive studies on neutrinos in cosmology can be found in Refs. [65−106].

Recently, the baryon acoustic oscillation (BAO) measurements from the second data release (DR2) of the Dark Energy Spectroscopic Instrument (DESI), derived from three years of observations, have been published. When combined with cosmic microwave background (CMB) and type Ia supernovae (SN) data, these measurements indicate a 2.8$ \sigma $–4.2$ \sigma $ preference for dynamical dark energy (DE) within the Chevallier–Polarski–Linder (CPL) model [107], providing stronger evidence than that obtained from the first data release (DR1) [108]. These significant deviations from the cosmological constant have drawn considerable attention within the cosmology community, motivating extensive studies on various aspects of cosmological physics, including the nature of DE and the properties of neutrinos [109−143]. In particular, Ref. [140] has recently reported that, in the $ \Lambda $CDM model, when BAO data from DESI DR1 are combined with SN data from DESY5 and CMB data from Planck and ACT, the cosmological upper limits on sterile neutrino parameters are constrained to $ N_{\rm{eff}} < 3.42 $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} < 0.217 $ eV at the 2$ \sigma $ level. In the CPL model, the corresponding constraints change to $ N_{\rm{eff}} < 3.37 $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} < 0.350 $ eV, indicating that the properties of DE can influence the cosmological constraints on sterile neutrino parameters.

In parallel, Ref. [39] conducted a preliminary investigation of sterile neutrinos within the $ \Lambda $CDM model, using joint observations of BNS mergers detected by 3G GW detectors and short $ \gamma $-ray burst (GRB) observed by a THESEUS-like mission. They found that the GW data provide a tight upper limit on $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $ and show a mild preference for $ \Delta N_{\rm{eff}} > 0 $ at the 1.1$ \sigma $ level, indicating that GW standard siren observations can significantly tighten the constraints on sterile neutrino parameters.

Taken together, these insights underscore the importance and timeliness of investigating sterile neutrinos within dynamical DE models using future GW observations. In this paper, we consider three representative dynamical DE models, namely the wCDM, holographic DE (HDE), and CPL models. We forecast the prospects for detecting sterile neutrinos in these models using joint GW-GRB observations. The primary aim of this study is to evaluate how future GW standard sirens can improve cosmological constraints on sterile neutrino parameters in dynamical DE models.

This work is organized as follows. In Sec. II, we introduce the models, cosmological data, and GW simulation used in this work. In Sec. III, we present the constraint results of cosmological parameters from the GW standard siren observations. The conclusion is given in Sec. IV.

In this section, we first introduce the cosmological models used in this study, followed by a brief overview of the cosmological datasets utilized, and finally describe the approach for simulating the GW data.

The baseline model considered is the $ \Lambda $CDM model. The extended model, $ \Lambda $CDM with sterile neutrinos ($ \Lambda $CDM+$ \nu_s $) includes eight base parameters, including the baryon density parameter $ \omega_b\equiv \Omega_b h^2 $, the cold dark matter density parameter $ \omega_c\equiv \Omega_c h^2 $, the ratio between the sound horizon and the angular diameter distance at decoupling $ \Theta_{\rm{s}} $, the optical depth to the reionization of the universe $ \tau $, the amplitude of initial curvature perturbation $ A_s $, the power-law spectral index $ n_s $, the effective number of relativistic species $ N_{\rm{eff}} $, and the effective sterile neutrino mass $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $. Note that the true mass of thermally distributed sterile neutrinos is given by $ m_{\rm{sterile}}^{\rm{thermal}}=(N_{\rm{eff}}-N_{\rm{SM}})^{-3/4}m_{\nu,{\rm{sterile}}}^{\rm{eff}} $, where $ N_{\rm{SM}} = 3.044 $ is the effective number of active neutrinos in the Standard Model framework [144−146]. In addition, the total mass of active neutrinos is fixed at $ \sum m_\nu = 0.06 $ eV in this study.

In our analysis, we consider three dynamical DE models, with their respective parameters detailed in Table 1, and followed by a brief introduction to each model.

Model	Parameter #	Parameters
wCDM+$ \nu_s $	9	$ \omega_b $, $ \omega_c $, $ \Theta_{\rm{s}} $, $ \tau $, $ \ln (10^{10}A_s) $, $ n_s $, $ N_{\rm{eff}} $, $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $, w
HDE+$ \nu_s $	9	$ \omega_b $, $ \omega_c $, $ \Theta_{\rm{s}} $, $ \tau $, $ \ln (10^{10}A_s) $, $ n_s $, $ N_{\rm{eff}} $, $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $, c
CPL+$ \nu_s $	10	$ \omega_b $, $ \omega_c $, $ \Theta_{\rm{s}} $, $ \tau $, $ \ln (10^{10}A_s) $, $ n_s $, $ N_{\rm{eff}} $, $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $, $ w_0 $, $ w_a $

Table 1.  Cosmological models and parameters used in this work.

● The wCDM+$ \nu_s $ model: Building on the $ \Lambda $CDM+$ \nu_s $ model, we consider the case where the equation-of-state (EoS) parameter w is treated as a constant free parameter.

● The HDE+$ \nu_s $ model: Extending the $ \Lambda $CDM+$ \nu_s $ model, we adopt the HDE model [147], which introduces an additional parameter c through the energy density definition $ \rho_{\rm{de}} = 3 c^2 M^2_{\rm{Pl}} L^{-2} $, where $ M_{\rm{Pl}} $ is the reduced Planck mass and L is the future event horizon of the universe. The parameter c is dimensionless and phenomenological, playing an important role in determining the properties of DE in this model. The EoS is given by $ w = - \dfrac{1}{3} - \dfrac{2}{3} \dfrac{\sqrt{\Omega_{\rm{de}}}}{c} $. For more details of the HDE model, see, e.g., Refs. [148−151].

● The CPL+$ \nu_s $ model: This model further extends $ \Lambda $CDM+$ \nu_s $ by introducing a time-varying EoS via the CPL parametrization [152, 153], $ w(a) = w_0 + w_a (1 - a) $, where $ w_0 $ denotes the present-day value of w and $ w_a $ characterizes its evolution with the scale factor a.

The cosmological observational data used in this study include the CMB, BAO, and SN. For the CMB data, we use the Planck TT, TE, EE spectra at $ \ell\geq 30 $, the low-$ \ell $ temperature commander likelihood, and the low-$ \ell $ SimAll EE likelihood from the Planck 2018 release [154]. For the BAO data, we adopt measurements from 6dFGS [155], SDSS-MGS [156], and BOSS-DR12 [157]. For the SN data, we employ the Pantheon sample, which consists of 1048 data points from the Pantheon complation [158]. For simplicity, we use "CBS" to denote the joint CMB+BAO+SN dataset.

Note that the primary objective of this work is to investigate the impact of joint GW and GRB observations on the cosmological constraints of sterile neutrino parameters within various dynamical DE models. To facilitate a direct and consistent comparison with previous studies on sterile neutrinos in the $ \Lambda $CDM model [39], we deliberately adopt the same combination of observational datasets (CBS). This approach allows us to isolate and evaluate the specific contribution of GW observations to the improvement in sterile neutrino constraints.

In this subsection, we present the method of simulating the joint GW standard sirens and GRB events. We focus on the synergy between the THESEUS-like GRB detector and the 3G GW observations. We utilize the simulation method as prescribed in Refs. [32, 36, 39, 40].

Based on the star formation rate [14, 159, 160], the merger rate of the BNS density with redshift in the observer frame is

where $ {\rm{d}}V(z)/{\rm{d}}z $ represents the comoving volume element, and $ {\cal{R}}_{\rm{m}}(z) $ is the BNS merger rate in the source frame, defined as

Here, $ t_{\rm{max}}=t_{\rm{H}} $ is the maximum delay time, $ t_{\rm{min}}=20 $ Myr is the minimum delay time, $ {\cal{R}}_{\rm{f}} $ is the cosmic star formation rate in the source frame based on the Madau-Dickinson model [161], $ t(z) $ is the age of the universe at the time of merger, $ t_{\rm{d}} $ is the delay time between BNS system formation and merger, and $ P(t_{\rm{d}}) $ represents the time delay distribution of $ t_{\rm{d}} $, for which we use the exponential time delay model [159], expressed as

with an e-fold time of $ \tau=100 $ Myr for $ t_{\rm{d}}>t_{\rm{min}} $.

In our analysis of BNS mergers, we adopt a local comoving merger rate of $ {\cal{R}}_{\rm{m}}(z=0)=920\; \rm Gpc^{-3}\; yr^{-1} $, which corresponds to the median estimate derived from the O1 LIGO and O2 LIGO/Virgo observation runs [162] and remains consistent with the findings from the O3 observation run [163]. We simulate a catalog of BNS mergers for a 10-year observation. For each source, the location $ (\theta,\phi) $, the polarization angle $ \psi $, the orientation angle $ \iota $, and the coalescence phase $ \psi_{\rm{c}} $ are all drawn from uniform distributions. To simplify, we adopt a Gaussian mass distribution. This distribution has a mean NS mass of $ 1.33\; M_{\odot} $ and a standard deviation of $ 0.09\; M_{\odot} $, where $ M_{\odot} $ denotes the solar mass [164, 165].

In the framework of the stationary phase approximation [166], the Fourier transform of the frequency-domain GW waveform for a detector network (with N detectors) is expressed as [167, 168]

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

where $ S_{\rm{n,k}}(f) $ is the one-side noise power spectral density of the k-th detector.

In this study, we consider the waveform in the inspiralling stage for the non-spinning BNS system. We employ the restricted Post-Newtonian approximation up to 3.5 PN order [169, 170]. The Fourier transform of the GW waveform for the k-th detector is expressed as

where the detailed forms of $ \Psi(f/2) $ and $ \varphi_{k,(2,0)} $ can be found in Refs. [168, 169], and the Fourier amplitude can be written as

where $ d_{\rm{L}} $ is the luminosity distance to the GW source, $ F_{+,k} $ and $ F_{\times,k} $ represent the antenna response functions of the k-th GW detector, $ {\cal{M}}_{\rm{chirp}}=(1+z)\eta^{3/5}M $ is the chirp mass of the binary system, $ \eta=m_1 m_2/M^2 $ is the symmetric mass ratio, and $ M=m_1+m_2 $ is the total mass of the binary system with $ m_1 $ and $ m_2 $ as the masses of the components. We use the GW waveform in the frequency domain, where time t is replaced by $ t_{f}=t_{\rm{c}}-(5 / 256) {\cal{M}}_{\rm{chirp}}^{-5 / 3}(\pi f)^{-8 / 3} $ [168, 169], with $ t_{\rm{c}} $ being the coalescence time.

Following the simulation of the GW catalog, we calculate the signal-to-noise ratio (SNR) for each GW event. Here, we adopt the SNR threshold to be 12 in our simulation. For low-mass systems, the combined SNR for the detection network of N independent detectors can be calculated by

The inner product is defined as

where $ {\boldsymbol{a}} $ and $ {\boldsymbol{b}} $ are column matrices with the same dimension, $ * $ denotes conjugate transpose, the lower cutoff frequency is set to $ f_{\rm{lower}}=1 $ Hz for ET and $ f_{\rm{lower}}=5 $ Hz for CE, and $ f_{\rm{upper}}=2/(6^{3/2}2\pi M_{\rm{obs}}) $ represents the frequency in the last stable orbit, with $ M_{\rm{obs}}=(m_1+m_2)(1+z) $.

For the short GRB model, we use the Gaussian structured jet profile model, which is based on the analysis of GW170817/GRB170817A [171],

where $ \theta_{\rm{v}} $ is the viewing angle, $ L_{\rm{iso}}(\theta_{\rm{v}}) $ represents the isotropically equivalent luminosity of short GRB observed at different values of $ \theta_{\rm{v}} $, $ L_{\rm{on}}=L_{\rm{iso}}(0) $ is the on-axis isotropic luminosity, $ \theta_{\rm{c}}=4.7^{\circ} $ is the characteristic core angle, and the direction of the jet is assumed to be aligned with the binary's orbital angular momentum, namely $ \iota=\theta_{\rm{v}} $.

For the distribution of the short GRB, we assume an empirical broken-power-law luminosity function, expressed as

where L is the peak luminosity in the 1–10000 keV energy range in the rest frame assuming isotropic emission, $ L_{*} $ is the characteristic parameter separating the two regimes, $ \alpha $ and $ \beta $ are the characteristic slopes that describe these regimes, respectively. Following Ref. [172], we use $ L_{*}=2\times10^{52} $ erg sec^-1, $ \alpha=-1.95 $, and $ \beta=-3 $. We adopt a standard low end cutoff in luminosity of $ L_{\rm{min}} = 10^{49} $ erg sec^-1, and we treat the on-axis isotropic luminosity $ L_{\rm{on}} $ as the peak luminosity L [14, 32, 160, 173]. For the THESEUS mission [174], a GRB detection is recorded if the value of observed flux is larger than the flux threshold $ P_{\rm{T}}=0.2$ ph s^-1 cm^-2 in the 50–300 keV band. For the GRB detection, we take a sky coverage fraction of 0.5 and a duty cycle of 80% [174]. Thus from the GW catalog that exceeds the threshold of 12, we can select the GW-GRB events based on the probability distribution $ \Phi(L){\rm{d}}L $.

For a network with N independent interferometers, the Fisher information matrix is defined as

where $ \theta_i $ represents nine GW parameters ($ d_{\rm{L}} $, $ {\cal{M}}_{\rm{chirp}} $, $ \eta $, $ \theta $, $ \phi $, $ \iota $, $ t_{\rm{c}} $, $ \psi_{\rm{c}} $, $ \psi $) for each GW event. The covariance matrix is equal to the inverse of the Fisher matrix, i.e., $ {\rm{Cov}}_{ij}=(F^{-1})_{ij} $. Therefore, the instrumental error $ \sigma_{d_{\rm{L}}}^{\rm{inst}} $ of the GW parameter $ \theta _i $ is given by $ \Delta\theta_i=\sqrt{{\rm{Cov}}_{ii}} $.

The total uncertainty of the luminosity distance $ d_{\rm{L}} $ can be written as

where the instrumental error of $ \sigma_{d_{\rm{L}}}^{\rm{inst}} $ is estimated by the Fisher information matrix, $ \sigma_{d_{\rm{L}}}^{\rm{lens}} $ is the weak-lensing error, and $ \sigma_{d_{\rm{L}}}^{\rm{pv}} $ is the peculiar velocity error.

The $ \sigma_{d_{\rm{L}}}^{\rm{lens}} $ is adopted from Refs. [175, 176],

The $ \sigma_{d_{\rm{L}}}^{\rm{pv}} $ of the GW source is given by [177]

where c is the speed of light in vacuum, $ H(z) $ is the Hubble parameter, and $ \sqrt{\langle v^2\rangle} $ is the peculiar velocity of the GW source, set to $ \sqrt{\langle v^2\rangle}=500\ {\rm km\ s^{-1}} $.

For the GW standard siren observation with N data points, the $ \chi^2 $ function is given by

where $ \vec{\Omega} $ is the set of cosmological parameters, $ {z}_i $, $ {d}_L^i $, and $ {\sigma}_{d_{\rm{L}}}^i $ are the i-th GW redshift, luminosity distance, and the total error of the luminosity distance, respectively.

In this work, we present the forecast for the search for sterile neutrinos in dynamical DE models using joint GW-GRB observations. We use the public Markov-chain Monte Carlo (MCMC) package $\mathtt{CosmoMC}$ [178] to infer the posterior probability distributions of the cosmological parameters. To demonstrate the impact of simulated GW data on constraining sterile neutrino parameters, we will analyze two scenarios of 3G GW observations: a single ET and the ET-CE-CE network, which consists of one ET detector and two CE-like detectors (one in the United States with a 40 km arm length and another one in Australia with a 20 km arm length, hereafter referred to as ET2CE). We utilize the ET sensitivity curve provided in Ref. [179] and the CE sensitivity curves from Ref. [180], as presented in Fig. 1 of Ref. [39]. For the GW detectors, due to the large uncertainty in the duty cycle, we consider only the ideal scenario in which all detectors are assumed to have a 100% duty cycle, as discussed in Ref. [181]. The geometry of the GW detector, characterized by parameters such as latitude ($ \varphi $), longitude ($ \lambda $), opening angle ($ \zeta $), and arm bisector angle ($ \gamma $), is described in detail in Table I of Ref. [39]. To ensure a direct comparison of constraints within the same parameter space, we adopt the best-fit cosmological parameter values from the CBS dataset as the fiducial values for simulating GW data in dynamical DE models. Thus, in our analysis, we use three data combinations: (1) CBS, (2) CBS+ET, and (3) CBS+ET2CE. We will report the constraint results in the next section.

Figure 1.  (color online) he triangular plots of the marginalized posterior distributions for cosmological parameters in the wCDM+$ \nu_s $ model, using the CBS, CBS+ET, and CBS+ET2CE data combinations, respectively.

Figure 2.  (color online) The triangular plots of the marginalized posterior distributions for cosmological parameters in the HDE+$ \nu_s $ model, using the CBS, CBS+ET, and CBS+ET2CE data combinations, respectively.

Figure 3.  (color online) The triangular plots of the marginalized posterior distributions for cosmological parameters in the CPL+$ \nu_s $ model, using the CBS, CBS+ET, and CBS+ET2CE data combinations, respectively.

Figure 4.  (color online) Two-dimensional marginalized posterior contours (1$ \sigma $ and 2$ \sigma $) in the $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $–$ N_{\rm{eff}} $ plane for the $ \Lambda $CDM+$ \nu_s $, wCDM+$ \nu_s $, HDE+$ \nu_s $, and CPL+$ \nu_s $ models using the CBS (a), CBS+ET (b), and CBS+ET2CE (c) data combinations, respectively.

In this section, we present the fitting results for the wCDM+$ \nu_s $, HDE+$ \nu_s $, and CPL+$ \nu_s $ models, respectively. We constrain these cosmological models using future GW standard siren observations. It is found that approximately 400 and 640 standard sirens can be detected over 10 years by the ET+THESEUS and ET2CE+THESEUS networks, respectively, with their redshift distributions provided in Ref. [39]. The fitting results are presented in Figs. 1–4 as well as Table 2. Note that we quote $ \pm 1\sigma $ errors for the parameters, but for the parameters $ N_{\rm{eff}} $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $, when central values cannot be determined, we instead quote the $ 2\sigma $ upper limits. We use $ \sigma(\xi) $ and $ \varepsilon(\xi)=\sigma(\xi)/\xi $ to represent the absolute error and the relative error of a parameter $ \xi $, respectively. For direct comparison, we also reproduce the constraint of the $ \Lambda $CDM+$ \nu_s $ [39] under the CBS, CBS+ET and CBS+ET2CE data combinations in Table 3.

	wCDM+$ \nu_s $				HDE+$ \nu_s $				CPL+$ \nu_s $
	CBS	CBS+ET	CBS+ET2CE		CBS	CBS+ET	CBS+ET2CE		CBS	CBS+ET	CBS+ET2CE
$ \sigma(\Omega_m) $	0.0078	0.0042	0.0042		0.0077	0.0040	0.0038		0.0081	0.0049	0.0048
$ \sigma(H_0) $	0.9850	0.0655	0.0655		1.1300	0.0750	0.0735		0.9700	0.0855	0.0880
$ \sigma(w/c/w_0) $a	0.0360	0.0160	0.0160		0.031	0.015	0.015		0.087	0.031	0.031
$ \sigma(w_a) $	...	...	...		...	...	...		0.39	0.18	0.18
$ \varepsilon(\Omega_m) $	2.542%	1.375%	1.375%		2.516%	1.309%	1.243%		2.620%	1.589%	1.556%
$ \varepsilon(H_0) $	1.432%	0.096%	0.096%		1.646%	0.109%	0.107%		1.409%	0.124%	0.128%
$ \varepsilon(w/c/w_0) $	3.472%	1.540%	1.540%		4.851%	2.344%	2.344%		9.315%	3.326%	3.330%
$ \varepsilon(w_a) $	...	...	...		...	...	...		78.000%	35.294%	35.294%
$ N_{\rm{eff}} $	<3.340	$ 3.142^{+0.032}_{-0.091} $	$ 3.143^{+0.035}_{-0.087} $		<3.562	$ 3.217^{+0.068}_{-0.062} $	3.219±0.058		<3.404	$ 3.165^{+0.045}_{-0.103} $	$ 3.162^{+0.044}_{-0.098} $
$ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $	<0.651	<0.565	<0.564		<0.479	<0.121	<0.115		<0.707	<0.691	<0.680
a w is for wCDM+$ \nu_s $, c is for HDE+$ \nu_s $ and $ w_0 $ is for CPL+$ \nu_s $.

Table 2.  Fitting results of the wCDM+$ \nu_s $, HDE+$ \nu_s $ and CPL+$ \nu_s $ models by using the CBS, CBS+ET, and CBS+ET2CE data combinations. We quote $ \pm 1\sigma $ errors for the parameters, but for the parameters $ N_{\rm{eff}} $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $, when central values cannot be determined, we quote the $ 2\sigma $ upper limits. Here, $ H_0 $ is in units of $ {\rm{km}}\ {\rm{s}}^{-1}\ {\rm{Mpc}}^{-1} $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $ is in units of eV. For a parameter $ \xi $, $ \sigma(\xi) $ and $ \varepsilon(\xi)=\sigma(\xi)/\xi $ represent its absolute and relative errors, respectively.

Parameter	CBS	CBS+ET	CBS+ET2CE
$ \sigma(\Omega_m) $	0.00640	0.00420	0.00405
$ \sigma(H_0) $	0.7500	0.0540	0.0520
$ \varepsilon(\Omega_m) $	2.055%	1.357%	1.308%
$ \varepsilon(H_0) $	1.101%	0.079%	0.076%
$ N_{\rm{eff}} $	<3.3446	$ 3.148^{+0.039}_{-0.094} $	$ 3.150^{+0.042}_{-0.089} $
$ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $	<0.5789	<0.4842	<0.4226

Table 3.  Same as Table 2, but in the $ \Lambda $CDM+$ \nu_s $ model [39].

First, we examine how the simulated GW data and the properties of DE affect the constraints on sterile neutrino parameters.

In Table 2, using CBS data, we obtain $ N_{\rm{eff}}<3.340 $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}}<0.651 $ eV for wCDM+$ \nu_s $, $ N_{\rm{eff}}<3.562 $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}}<0.479 $ eV for HDE+$ \nu_s $, $ N_{\rm{eff}}<3.404 $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}}<0.707 $ eV for CPL+$ \nu_s $. Using CBS+ET data, we obtain $ N_{\rm{eff}}=3.142^{+0.032}_{-0.091} $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}}<0.565 $ eV for wCDM+$ \nu_s $, $ N_{\rm{eff}}=3.217^{+0.068}_{-0.062} $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}}<0.121 $ eV for HDE+$ \nu_s $, $ N_{\rm{eff}}=3.165^{+0.045}_{-0.103} $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}}<0.691 $ eV for CPL+$ \nu_s $. Using CBS+ET2CE data, we obtain $ N_{\rm{eff}}=3.143^{+0.035}_{-0.087} $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}}<0.564 $ eV for wCDM+$ \nu_s $, $ N_{\rm{eff}}=3.219\pm0.058 $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}}<0.115 $ eV for HDE+$ \nu_s $, $ N_{\rm{eff}}=3.162^{+0.044}_{-0.098} $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}}<0.680 $ eV for CPL+$ \nu_s $.

We find that the CBS dataset alone does not enable a precise determination of either parameter $ N_{\rm{eff}} $ or $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $, regardless of the dynamical DE model considered; only upper limits can be obtained for both parameters. However, when GW standard sirens are incorporated into the data combination, the constraints on $ N_{\rm{eff}} $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $ are significantly improved across all dynamical DE models.

The preference for $ \Delta N_{\rm{eff}}>0 $ reaches approximately $ 1.1\sigma $ (CBS+ET and CBS+ET2CE) for the wCDM+$ \nu_s $ model, $ 2.8\sigma $ (CBS+ET) and $ 3.0\sigma $ (CBS+ET2CE) for the HDE+$ \nu_s $ model, and $ 1.2\sigma $ (CBS+ET and CBS+ET2CE) for the CPL+$ \nu_s $ model. Moreover, the inclusion of GW data helps reduce the $ 2\sigma $ upper limit on $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $ by 13.2%–13.4%, 74.7%–76.0%, and 2.3%–3.8% for the wCDM+$ \nu_s $, HDE+$ \nu_s $, and CPL+$ \nu_s $ models, respectively. These results clearly demonstrate that GW observations can significantly tighten the constraints on both $ N_{\rm{eff}} $ and $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $ (see also Figs. 1–3). This finding is consistent with previous studies on sterile neutrinos in the $ \Lambda $CDM model using GW data [39].

To explicitly illustrate how the properties of DE influence the constraints on sterile neutrino parameters, Fig. 4 shows the joint constraints on the $ \Lambda $CDM+$ \nu_s $, wCDM+$ \nu_s $, HDE+$ \nu_s $, and CPL+$ \nu_s $ models in the $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $–$ N_{\rm{eff}} $ plane, derived from the CBS, CBS+ET, and CBS+ET2CE data combinations. We find that, compared to the $ \Lambda $CDM+$ \nu_s $ model (see also Table 3), in the wCDM+$ \nu_s $ and CPL+$ \nu_s $ models the limits on $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $ become much looser, but in the HDE+$ \nu_s $ model the limits become much more stringent. This clearly indicates that the DE properties can evidently affect the upper limits on $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $.

Therefore, both GW data and the properties of DE play an important role in constraining the sterile neutrino parameters, with GW observations significantly enhancing the precision of these constraints.

Next, we discuss how the GW data help to improve the constraint accuracies for other cosmological parameters.

For the wCDM+$ \nu_s $ model, Fig. 1 presents the constraint results from the CBS, CBS+ET, and CBS+ET2CE combinations. It is evident that incorporating GW data significantly improves the constraints on $ \Omega_m $, $ H_0 $, and w. Specifically, compared to CBS alone, the inclusion of ET or ET2CE data improves the constraints on $ \Omega_m $, $ H_0 $, and w by 46.2%, 93.4%, and 55.6%, respectively.

For the HDE+$ \nu_s $ model, Fig. 2 presents the corresponding constraint results. The inclusion of GW data also leads to notable improvements in the constraints on $ \Omega_m $, $ H_0 $, and the model parameter c. Specifically, compared to CBS alone, the inclusion of ET data improves the constraints on $ \Omega_m $, $ H_0 $, and c by 48.1%, 93.4% and 51.6%, respectively. When ET2CE data are included, the corresponding improvements reach 50.6%, 93.5% and 51.6%, respectively.

For the CPL+$ \nu_s $ model, Fig. 3 shows the constraints from the CBS, CBS+ET, and CBS+ET2CE combinations. It is clear that incorporating GW data significantly improves the constraints on $ \Omega_m $, $ H_0 $, $ w_0 $, and $ w_a $. With the inclusion of ET data, the constraints on $ \Omega_m $, $ H_0 $, $ w_0 $, and $ w_a $ are improved by 39.5%, 91.2%, 64.4%, and 53.8%, respectively, and by 40.7%, 90.9%, 64.4%, and 53.8% with the ET2CE data.

These results demonstrate that incorporating GW data significantly enhances the precision of cosmological parameter constraints, including $ \Omega_m $, $ H_0 $, w, c, $ w_0 $, and $ w_a $.

This work aims to forecast the potential for detecting sterile neutrinos in dynamical DE cosmologies using joint GW-GRB observations. We consider three representative models, the wCDM, HDE, and CPL models, to investigate how future GW data may improve constraints on sterile neutrino and other cosmological parameters. Two GW observation strategies are explored: the ET alone and the ET2CE network. To evaluate the impact of GW data, we also incorporate existing CMB, BAO, and SN data for comparative and joint analyses.

We find that the properties of DE could significantly influence the constraints on the sterile neutrino parameters. Compared to the $ \Lambda $CDM model, in the HDE model the limits on $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $ become much more stringent, but in the wCDM model and the CPL model the limits become much looser.

In addition, we find that the inclusion of GW data significantly improves the constraints on sterile neutrino parameters in dynamical DE models, compared to the current limits derived from CMB+BAO+SN observations. Specifically, the GW data significantly tighten the constraint on $ N_{\rm{eff}} $ and indicate a preference for $ \Delta N_{\rm{eff}}>0 $ at approximately 1.1$ \sigma $ in the wCDM model, 3.0$ \sigma $ in the HDE model, and 1.2$ \sigma $ in the CPL model. Moreover, the upper limits on $ m_{\nu,{\rm{sterile}}}^{\rm{eff}} $ are reduced by 13.2%–13.4% for wCDM, 74.7%–76.0% for HDE, and 2.3%–3.8% for CPL, respectively. These results demonstrate that GW observations can substantially enhance the prospects for detecting sterile neutrinos.

Furthermore, we find that the GW data also play a significant role in improving the constraints on other cosmological parameters, including $ \Omega_m $, $ H_0 $, and the EoS parameters of DE, namely, w in the wCDM model, c in the HDE model, and $ w_0 $ and $ w_a $ in the CPL model.

These findings highlight that GW observations significantly improve the constraints on most cosmological parameters and enhance the prospects for detecting sterile neutrinos within various DE models. With upgrades to GW detectors, a substantial sample of standard sirens is expected to be produced, which will play a crucial role in further advancing sterile neutrino searches.

We thank Tian-Nuo Li and Guo-Hong Du for their helpful discussions.