Prospects for Probing the Interaction between Dark Energy and Dark Matter Using Gravitational-wave Dark Sirens with Neutron Star Tidal Deformation

Gravitational wave (GW) standard siren observations provide a rather useful tool to explore the evolution of the Universe. In this work, we wish to investigate whether dark sirens with neutron star (NS) deformation from third-generation GW detectors could help probe the interaction between dark energy and dark matter. We simulate the GW dark sirens of four detection strategies based on 3 yr observation and consider four phenomenological interacting dark energy (IDE) models to perform cosmological analysis. We find that GW dark sirens could provide tight constraints on Ωm and H 0 in the four IDE models but do not perform well in constraining the dimensionless coupling parameter β in models of the interaction proportional to the energy density of cold dark matter. Nevertheless, the parameter degeneracy orientations of cosmic microwave background (CMB) and GW are almost orthogonal, and thus, the combination of them could effectively break cosmological parameter degeneracies, with the constraint errors of β being 0.00068–0.018. In addition, we choose three typical equations of state (EoSs) of an NS, i.e., SLy, MPA1, and MS1, to investigate the effect of an NS’s EoS on cosmological analysis. The stiffer EoS could give tighter constraints than the softer EoS. Nonetheless, the combination of CMB and GW dark sirens (using different EoSs of an NS) shows basically the same constraint results of cosmological parameters. We conclude that the dark sirens from 3G GW detectors would play a crucial role in helping probe the interaction between dark energy and dark matter, and the CMB+GW results are basically not affected by the EoS of an NS.

Among various extensions to the ΛCDM model, there is a class of models known as the interaction dark energy (IDE) models, in which the direct and nongravitational interaction between dark energy and dark matter is considered.On the one hand, the IDE models can not only help alleviate the Hubble tension (Guo et al. 2018;Feng et al. 2020b) but also alleviate the coincidence problem of dark energy (Cai & Wang 2005;Zhang 2005).On the other hand, the investigation of the IDE models can help probe the fundamental nature of dark energy and dark matter.Although CMB can give rather tight constraints on the ΛCDM model, the newly extended parameters beyond the ΛCDM model will severely degenerate with other parameters when the CMB data alone are used to constrain the extended cosmological models.Hence, other cosmological probes need to be combined with the CMB data to break the cosmological parameter degeneracies.
In the next decades, the third-generation (3G) detectors, the Einstein Telescope (ET; Punturo et al. 2010) and the Cosmic Explorer (CE; Abbott et al. 2017d) will be operational, with sensitivities improved one order of magnitude over the current GW detectors and a large number of GW events can be observed (Evans et al. 2021).Nevertheless, the detection of EM counterparts is rather difficult.Recent forecasts show that only ∼ 0.1% of binary neutron star (BNS) mergers detected by the 3G detectors can have detectable EM counterparts (Yu et al. 2021).Moreover, bright sirens can only give a tight constraint on the Hubble constant but do not perform well in constraining other cosmological parameters.
In general, only BNS and black hole (BH)-NS merger events have the potential for EM counterparts, and in most cases, these events are not observable in EM waves.Binary black hole (BBH) merger events do not have EM counterparts (it should be noted that the events of supermassive or massive BH binaries may have observable EM counterparts; see, e.g., Tamanini et al. 2016, Wang et al. 2020b, 2022a,b, and Zhao et al. 2020).For GW events without observable EM counterparts, there are actually methods to conduct cosmological studies using standard sirens.The commonly employed approach is to statistically determine the redshifts of GW events by cross-correlating them with galaxy catalogs; see, e.g., Fishbach et al. (2019), Soares-Santos et al. (2019), Palmese et al. (2020), Abbott et al. (2021aAbbott et al. ( , 2023a)), and Palmese et al. (2023).This type of GW standard siren is referred to as a dark siren.
In fact, it is possible to establish the distance-redshift relation solely through GW observations without the need for EM observations, including galaxy catalogs.Messenger & Read (2012) proposed that in BNS merger events, the additional tidal phase caused by the tidal deformation of an NS can aid in determining the redshifts of GW events.Consequently, both luminosity distance and redshift can be independently determined from GW observations alone.Since this method does not involve EM counterpart observations, in this paper, we also refer to such standard sirens as dark sirens.It is important to note that the dark sirens discussed in this paper are distinct from the dark sirens using the galaxy catalog method.
Therefore, employing this method, all the detected BNS mergers can be used to perform cosmological analysis; see, e.g., Del Pozzo et al. (2017), Wang et al. (2020a), Chatterjee et al. (2021), Dhani et al. (2022) Ghosh et al. (2022), andJin et al. (2023a) for recent works.However, the tidal deformation is related to the equation of state (EoS) of NS, which may lead to different constraint results.By solely using the dark sirens, the EoS of dark energy and the Hubble constant can be precisely measured, provided that the EoS of an NS is perfectly known (Jin et al. 2023a).
These dark sirens show a strong ability to constrain cosmological parameters.Therefore, in this work, we wish to investigate the potential of the dark sirens in constraining the IDE models.The aim of this work is to answer three key questions: (i) what precisions the cosmological parameters in the IDE models could be measured by solely using the dark sirens, (ii) how the GW dark sirens break the cosmological parameter degeneracies generated by the EM observations, and (iii) how the EoS of an NS affects the constraints on the cosmological parameters.Through this work, we wish to shed light on these aspects and explore the potential role of the dark sirens in probing the interaction between dark energy and dark matter.
This work is organized as follows.In Sec. 2, we introduce the methodology used in this work.In Sec. 3, we report the constraint results and make some relevant discussions.The conclusion is given in Sec. 4.

Interacting Dark Energy Models
In a flat Friedmann-Roberston-Walker universe, the Friedmann equation is given by where 3 2 pl  2 is the critical density of the universe, and  de ,  c ,  b , and  r represent the energy densities of dark energy, cold dark matter, baryon, and radiation, respectively.
In the IDE models, if assuming a direct interaction between dark energy and cold dark matter, the energy conservation equations can be given by where the dot is the derivative with respect to the cosmic time ,  is the Hubble parameter,  is the EoS of dark energy, and  is the energy transfer rate.In this work, we consider four typical phenomenological forms of , i.e.,  =   de (IΛCDM1),  =   c (IΛCDM2),  =  0  de (IΛCDM3), and  =  0  c (IΛCDM4).Here  is a dimensionless coupling parameter that describes the interaction strength between dark energy and dark matter. > 0 indicates cold dark matter decaying into dark energy,  < 0 indicates dark energy decaying into cold dark matter, and  = 0 indicates no interaction between dark energy and cold dark matter.Note that we wish to investigate the potential of GW dark sirens in probing the interaction (or constraining ), and thus we consider  = −1 for simplicity.
In the IDE models, there is a problem of early-time perturbation instability (Majerotto et al. 2009;Clemson et al. 2012) because the cosmological perturbations of dark energy in the IDE models will be divergent in a part of the parameter space, which leads to the ruin of the IDE cosmology in the perturbation level.To avoid this problem, Li et al. (2014Li et al. ( , 2016)); Li & Zhang (2023) extended the parameterized post-Friedmann (PPF) approach (Fang et al. 2008;Hu 2008) to the IDE models, referred to as the ePPF approach.This approach can

Note.
Here  is the latitude,  is the longitude,  is the orientation angle of the detector's arms measured counterclockwise from east of the Earth to the bisector of the interferometer arms, and  is the opening angle between the detector's two arms.
safely calculate the cosmological perturbations in the whole parameter space of the IDE models.In this work, we employ the ePPF approach to treat the cosmological perturbations; see e.g., Zhang (2017); Feng et al. (2018) for the application of the ePPF approach.

Mock Gravitational-wave Dark Siren Data
Under the stationary phase approximation (Zhang et al. 2017), the Fourier transform of the time-domain waveform for a detector network is given by Wen & Chen (2010); Zhao & Wen (2018) where Φ is the  ×  diagonal matrix with Φ   = 2     (n • r  ), n is the propagation direction of a GW, r  is the location of the th detector, and  is given by ĥ In this work, we consider the waveform in the inspiralling stage for the nonspinning BNS system.We adopt the restricted post-Newtonian (PN) approximation and calculate the waveform to the 3.5 PN order (Cutler et al. 1993;Blanchet & Iyer 2005;Sathyaprakash & Schutz 2009).The Fourier transform of the GW waveform of the th detector is given by Sathyaprakash & Schutz (2009) where the Fourier amplitude A  is given by    2009) where  L is the luminosity distance to the GW source,  +, and  ×, are the antenna response functions1 of the th GW detector which are related to the locations of the GW source and the GW detector (the adopted coordinates of GW observatories are listed in Table 1 and shown in Figure 1),  is the inclination angle between the binary's orbital angular momentum and the line of sight, M chirp = (1 + ) 3/5  is the observed chirp mass,  =  1 +  2 is the total mass of binary system with component masses  1 and  2 ,  =  1  2 /( 1 +  2 ) 2 is the symmetric mass ratio,  c is the coalescence phase, and the detailed forms of coefficients   could be found in e.g., Sathyaprakash & Schutz (2009).Note that the GW waveform in the frequency domain is adopted, i.e., time  is replaced by Maggiore ( 2007) where  c is the coalescence time.
The tidal contribution to the GW phase in the frequency domain is given by Messenger & Read (2012) where the sum is comprised of the components of the binary system,   =   /( 1 +  2 ), and  is the tidal deformation parameter, which is the function of the mass of an NS.Given that various EoSs of an NS can yield distinct redshift measurements, a stiffer EoS is deemed to provide more accurate redshift determinations (Messenger & Read 2012).In order to investigate the impact of an NS's EoS on the cosmological parameter estimations, we chose the SLy (Douchin & Haensel 2001), MPA1 (Müther et al. 1987), and MS1 (Mueller & Serot 1996)  For the BNS merger rate, we used the methodology detailed in Jin et al. (2023a).We assume the Madau-Dickinson star formation rate (Madau & Dickinson 2014), with an exponential time delay (an e-fold time of 100 Myr; Vitale et al. 2019).Meanwhile, we adopt the local comoving merger rate to be 320 Gpc −3 yr −1 , which is the estimated median rate based on the O3a observation (Abbott et al. 2021b).In Figure 2, we show the expected detection rates as a function of redshift for ET, CE, the 2CE network, and the ET2CE network.Note that in the following discussions, we consider a 3 yr observation of the four detection strategies.
Throughout this work, we simulate the BNS mergers with random binary orientations and sky directions.The sky direction, binary inclination, polarization angle, and coalescence phase are randomly chosen in the ranges of cos  ∈ For the mass of an NS, we adopt the mass distribution given by the latest LIGO-VIRGO-KAGRA observation (Abbott et al. 2023b).
After simulating the GW catalog, we need to calculate the detection number of GW events.The signal-to-noise ratio (S/N) is given by The inner product is defined as where * represents complex conjugate,  lower is the lower cutoff frequency,  upper = 2/ 6 3/2 2 obs is the frequency at the last stable orbit with  obs = ( 1 +  2 ) (1 + ), and  n (  ) is the one-side power spectral density (PSD).The adopted PSD forms for ET 2 and CE 3 are shown in Figure 3.Note that we take into account the triangular configuration of ET in the calculation of S/N.We adopt the threshold of S/N to be 8 in the simulation.
Then we use the Fisher information matrix to estimate the instrumental error of the luminosity distance.The Fisher information matrix of a network including  independent GW detectors is given by where  denotes 12 parameters ( L , , M chirp , ,  c ,  c , , , , , , ) for a given BNS system.The covariance matrix is equal to the inverse of the Fisher information matrix, and thus we can use the Fisher information matrix to calculate the measurement errors of the above 12 GW parameters.The error of the parameter   can be calculated by where  is the total Fisher information matrix.
For the measurement errors of  L , we consider the instrumental error  inst  L and weak-lensing error  lens  L .The total error of  L is The instrumental error of  L could be obtained by Eq. ( 15), i.e.,  inst d L = Δ L = √︁ ( −1 ) 11 .The error caused by weak lensing is adopted from Hirata et al. (2010), Tamanini et al. (2016), Speri et al. (2021), For the measurement error of redshift Δ using the dark siren method, i.e., Δ = √︁ ( −1 ) 22 , we show the relative errors of redshift distributions for ET, CE, 2CE, and ET2CE in the left panel of Figure 4. We see that the measurement precisions of redshifts are mainly 10%-30%, which is consistent with the result given in Jin et al. (2023a)

Cosmological Data
For comparison, we also employ the current cosmological data, i.e., the CMB data, the baryon acoustic oscillation (BAO) data, and the type Ia supernova (SN) data.For the  (Brout et al. 2022).
We adopt the Markov Chain Monte Carlo method (Lewis & Bridle 2002) to maximize the likelihood and infer the posterior distributions of cosmological parameters.The detailed calculation of the likelihood can be referred to in Jin et al. (2023a).

RESULTS AND DISCUSSIONS
In this section, we shall report the constraint results of the cosmological parameters.We consider the IΛCDM1 ( =   de ), IΛCDM2 ( =   c ), IΛCDM3 ( =  0  de ), and IΛCDM4 ( =  0  c ) models to perform cosmological analysis.For the GW observations, we consider four detection strategies, i.e., ET, CE, the 2CE network, and the ET2CE network, to constrain the cosmological models.In order to show the ability of GW dark sirens to break cosmological parameter degeneracies, we also show the constraint results of CMB, CMB+ET, CMB+CE, CMB+2CE, CMB+ET2CE, and CMB+BAO+SN.The posterior distribution contours (1 and 2) for the interested cosmological parameters are shown in Figures 5-7 and summarized in Table 2.The evolutions of /( 0  cr0 ) in the IΛCDM1 and IΛCDM2 models are shown in Figure 8.In order to show the impact of different EoSs of an NS on estimating cosmological parameters using these dark sirens, we choose three typical EoSs of an NS to perform cosmological analysis and show the constraint results in Figure 9.

Constraints on the Interacting Dark Energy models
In Figure 5, we show the constraint results of the ET2CE network in the - 0 (left panel) and -Ω m (right panel) planes for the IΛCDM1, IΛCDM2, IΛCDM3, and IΛCDM4 models.From the left panel, we can see that IΛCDM2 gives the best constraints on cosmological parameters, followed by IΛCDM4, while IΛCDM1 and IΛCDM3 give comparable constraint results.Concretely, in the case of IΛCDM2, ET2CE gives () = 0.0130 and ( 0 ) = 0.071 km s −1 Mpc −1 , which is 43.5% and 11.3% better than those of IΛCDM1.From the right panel, we can see that the above results still hold.As shown in Table 2, while ET2CE gives the best constraint results among ET, CE, 2CE, and ET2CE, it still does not performs well in constraining the dimensionless coupling parameter , whereas GW dark sirens can all precisely measure the Hubble constant.Since IΛCDM2 gives the best constraint results among the IDE models, we choose IΛCDM2 as a representative to make the following discussion.
In Figure 6, in order to show the ability of GW dark sirens to break the cosmological parameter degeneracies, we show the constraints in the - 0 and -Ω m planes for the IΛCDM2 model.We can see that the parameter degeneracy orientations of CMB and ET2CE are almost orthogonal, and thus, their combination can effectively break the cosmological parameter degeneracies.Furthermore, the contours of CMB+ET2CE are much smaller than those of CMB+BAO+SN, indicating that the ability of dark sirens from ET2CE to break the cosmological parameters is much better than that of the current BAO and SN data.Concretely, CMB+ET2CE gives () = 0.00068, ( 0 ) = 0.041 km s −1 Mpc −1 , and (Ω m ) = 0.00087, which are 70.4%,97.6%, and 96.5% better than those of CMB, respectively.Moreover, the constraints on ,  0 , and Ω m using CMB+ET2CE are 43.3%, 93.7%, and 89.3% are better than those using CMB+BAO+SN, respectively.
In Figure 7, we show the constraint results of CMB+ET, CMB+CE, CMB+2CE, and CMB+ET2CE for the IΛCDM2 model.We can see that CMB+ET2CE also gives the best constraint results, followed by CMB+2CE, CMB+CE, and CMB+ET, but the difference is not as big as the constraint results of the single GW dark sirens, as shown in Table 2.Moreover, we can see that even CMB+ET gives much better constraint results than the current CMB+BAO+SN data, showing the potential of the GW dark sirens in breaking the cosmological parameter degeneracies.
In Figure 8, we show the evolutions of /( 0  cr0 ) versus  in the IΛCDM1 and IΛCDM2 models.Here  cr0 is the present-day critical density of the universe.

Impact of the Equation of State of Neutron Stars on Performing Cosmological Analysis
In this subsection, we wish to investigate the impact of an NS's EoS on estimating cosmological parameters.As pointed out in Messenger & Read (2012)  an NS could lead to different measurement precision of redshift when measuring the redshift from the tidal deformation of an NS.Hence, we wish to offer an answer to the question of what impact an NS's EoS will have on estimating the cosmological parameters.Here we choose three typical (soft, intermediate, and stiff) EoSs of an NS, i.e., SLy, MPA1, and MS1, to perform cosmological analysis.Note that the following discussions are based on the IΛCDM2 model using ET2CE.See Figure 9 and Table 3 for the results of the cosmologicalconstraints.
In the left panel of Figure 9, we show the constraint results of ET2CE based on SLy, MPA1, and MS1.We find that the stiffer MS1 gives tight constraints than the softer MPA1 and SLy.This is because the stiffer EoS can lead to more precise redshift measurement, as shown in the right panel of Figure 4, which is also consistent with the redshift measurement forecasted in Messenger & Read (2012).As mentioned above, the GW dark sirens can effectively break the cosmological parameter degeneracies generated by the EM observation.In the right panel of Figure 9, we show the constraint results of the combination of CMB and GW dark sirens.We find that the combination of CMB and GW dark sirens shows basically the same constraint results.This means that in the IDE scenario, an NS's EoS does have an impact on estimating cosmological parameters, whereas the combination results of  CMB and GW dark siren are basically not affected by the EoS of an NS.

CONCLUSION
In this work, we investigate the potential of the GW dark sirens from the 3G GW detectors in probing the interaction between dark energy and dark matter.Based on the 3 yr observation, we consider four detection strategies, i.e., ET, CE, the 2CE network, and the ET2CE network, to perform cosmological analysis.Four phenomenological IDE models are considered, i.e., IΛCDM1 ( =   de ), IΛCDM2 ( =   c ), IΛCDM3 ( =  0  de ), and IΛCDM4 ( =  0  c ).
We find that GW dark sirens can provide tight constraints on Ω m and  0 in all the IDE models but do not perform well in constraining  in the IΛCDM2 model.In the other IDE models, the constraints on  are all better than the result of the current CMB+BAO+SN data.Nonetheless, we find that in the IΛCDM2 model, the parameter degeneracy orientations of CMB and GW dark sirens are almost orthogonal, and thus, the combination of them could effectively break the cosmological parameter degeneracies and significantly improve the measurement precisions of the cosmological parameters.With the addition of the GW dark sirens to the CMB data, the constraints on cosmological parameters could be improved by 70.4%-97.6%.CMB+ET2CE gives () = 0.00068 in the IΛCDM2 model, which is much better than that of the current CMB+BAO+SN data.
In addition, we investigate the impact of the EoS of an NS on estimating cosmological parameters.The stiffer EoS could give tighter constraints than the softer EoS.However, the combination of CMB and GW dark sirens based on different EoSs of an NS shows basically the same constraint results, implying that the EoS of an NS has almost no impact on estimating cosmological parameters in the case of the combination of CMB and GW dark sirens.Our results show that the GW dark sirens from the 3G GW detectors are expected to play an important role in helping probe the interaction between dark energy and dark matter.

Figure 1 .Figure 2 .
Figure 1.Locations and orientations of the ground-based GW observatories.The second-generation GW detectors are shown in black lines (LIGO-Virgo-KAGRA Scientific Collaboration 2018; Ashton et al. 2019) and the 3G GW detectors considered in this work are shown in red lines (Di Giovanni et al. 2021; Borhanian 2021).

Figure 3 .
Figure 3.The sensitivity curves of ET and CE used in this work.Note that we consider one 20-km CE in Australia and one 40-km CE in the US.

Figure 4 .
Figure 4. Relative errors of redshift distributions.Left panel: distributions of Δ/ for ET, CE, 2CE, and ET2CE.Right panel: distributions of Δ/ for SLy, MPA1, and MS1 based on the observation of ET2CE.
models as representative EoSs of an NS, representing the soft, intermediate, and stiff EoSs of an NS, respectively.We use the linear function to fit the  −  relation in the mass range of [1.2, 2]  ⊙ , which is consistent with the current observations.We obtain  =  +  with  = −2.03and  = 4.43,  = −1.37 and  = 4.67, and  = −3.67 and  = 12.69 for the fitting of SLy, MPA1, and MS1.

Figure 8 .
Figure 8.The evolutions of /( 0  cr0 ) in the IΛCDM1 (left panel) and IΛCDM2 (right panel) models.The blue and pink shaded regions represent the 1 constraints from the CMB+BAO+SN and CMB+ET2CE data, respectively.The black dashed lines denote the ΛCDM model ( = 0).

Figure 9 .
Figure 9. Constraints on cosmological parameters using dark sirens from ET2CE based on the SLy, MPA1, and MS1 models.Left panel: constraints on the cosmological parameters using the dark sirens from ET2CE based on the SLy, MPA1, and MS1 models in the IΛCDM2 model.Right panel: constraints on the cosmological parameters using the CMB+ET2CE data based on the SLy, MPA1, and MS1 models in the IΛCDM2 model.

Table 2 .
The Absolute Errors (1) and the Relative Errors of the Cosmological Parameters.