Forecasts on interacting dark energy with standard sirens

We present the predictions with standard sirens at Gravitational Waves detectors, such as the Laser Interferometer Space Antenna (LISA) and the Einstein Telescope (ET), for interacting dark energy theories. We focus on four models characterised by couplings between the dark energy field and the dark matter fluid arising from conformal or disformal transformations of the metric, along with an exponential self-interacting potential. To this purpose we construct mock catalogues and perform a Markov Chain Monte Carlo analysis by considering ET and LISA standard sirens, and also their combination with Baryon Acoustic Oscillations (BAO) and Supernovae Ia (SNIa) data. We find that in all the four models considered, the accuracy on the $H_0$ parameter increases by one order of magnitude at 1$\sigma$ when compared to the SNIa+BAO data set, possibly shedding light in the future on the origin of the $H_0$-tension. The combination of standard sirens with SNIa+BAO allows to improve the accuracy on some coupling and exponential parameters, hinting at future prospects for constraining interactions in the dark sector.


I. INTRODUCTION
Understanding the nature of the dark sector of the Universe is one of the greatest endeavours of Cosmology at the present.This comprises the weakly interacting dark matter (DM) -responsible for the formation and dynamics of structures in the Universe -and dark energy (DE) -the driver of the late time cosmic acceleration.Together these components dominate about 95% of the energy budget of the Universe.In the Λ-cold-dark-matter (ΛCDM) scenario, i.e. the standard model of cosmology, DE is portrayed simply as a cosmological constant, Λ.This model comes with some theoretical issues [1][2][3], among which is the fact that fundamental theories do not properly account for the currently measured small value of the cosmological constant.ΛCDM also requires a primordial inflationary period to explain the geometrical flatness, Cosmic Microwave Background (CMB) smoothness, and initial conditions for large-scale structures.More recently, observational tensions on the value of the cosmological parameters H 0 [4][5][6][7][8][9][10] and σ 8 [11][12][13][14][15][16], measured by early-and late-Universe probes, increased the motivation to investigate alternative models of gravity [17].
Consequently, alternative theories are explored by cosmologists in which Λ is promoted to a dynamical DE scalar field, ϕ, namely quintessence [18,19] (see [20] for a review), which evolves in time according to its self-interaction potential.While in the standard ΛCDM scenario the two dark components do not directly couple with each other, in a dynamical DE model one can instead consider that they experience some non-minimal interaction.Such constructed models are referred to as coupled quintessence models [21,22].The dynamics of the field, along with the dark interaction, could provide a more natural explanation of the accelerated expansion, while also addressing the observational tensions [23].Nevertheless, the coupling can be formalised at the Lagrangian level through what is known as a conformal/disformal transformation of the metric tensor [24][25][26][27][28][29][30].If this transformation depends directly on the quintessence field, then this is physically equivalent to considering that the DM particles propagate on the geodesics of the transformed metric, ḡµν .In the conformal case, this is achieved from a re-scaling of the metric and, consequently, of time-and space-like norms and intervals alike, while preserving the light cones: where C is the conformal function.These find important applications in modified gravity theories as they preserve the structure of Scalar-Tensor theories of the Brans-Dicke form [31].Alternatively, one can consider that the metric transformation should depend on the first-order partial derivatives of the scalar field as well.This results in a disformal transformation: where C and D are the conformal and disformal functions, respectively and X = −g µν ∂ µ ϕ∂ ν ϕ/2.This gives rise to a more intricate scenario, with a distortion of the metric defined directionally according to the gradient of ϕ.First introduced by Bekenstein [25], the disformal transformations re-surged in the cosmological literature [27] when, in analogy to the conformal case, it was shown that they preserve a more general class of Scalar-Tensor theories categorised in the Horndeski Lagrangian [28].Disformal transformations in cosmology arise naturally in brane-world models [32,33] and have been the focus of many theoretical proposals for the nature of the dark sector and their interactions [26,[34][35][36][37][38][39][40][41][42].The stability conditions for the functions C and D have been discussed in Refs.[26][27][28]40] and the case in which D ≡ D(ϕ, X) has been discussed in Refs.[43,44].
The coupling in the dark sector gives rise to an additional gravitational fifth force in the Universe between DM particles, mediated by the DE field.This new force leaves distinct features in the background equations, as well as signatures in the cosmological density perturbations that describe the formation of structures [45].Although these deviations from the benchmark model are constrained to be small (especially at the background level), they are still expected to leave detectable, characteristic observational imprints that the data can probe.These are essential to test the viability of such alternative models by identifying the range of validity of the parameter space and the robustness of its predictions.
In the past few years, we have witnessed the rise of gravitational wave (GW) astronomy as a new independent probe of gravitational effects [46].An accurate redshift-luminosity relation can be constructed when GW events are combined with an electromagnetic (EM) counterpart multi-messenger signal.These observations become standard sirens [47], analogous to the standard candles used in local EM measurements.So far, only one GW event, GW170817, with a corresponding EM counterpart, GRB170817A, has been detected by the Laser Interferometer Gravitational-Wave Observatory (LIGO)-Virgo and the International Gamma-ray Astrophysics Laboratory (INTEGRAL)-Fermi collaborations, respectively, and which originated from the merger of a binary pair of neutron stars [48,49].This single combined detection had a strong impact on the allowed modifications to the gravitational interaction by ruling out many proposals [50][51][52][53][54] with many other models further constrained [55][56][57][58][59].
Current GW detectors, (advanced) Virgo [60], (advanced) LIGO [61] and the Kamioka Gravitational Wave Detector (KAGRA) [62], are secondgeneration (2G) ground-based detectors, with another one under planning (2030), the Indian Initiative in Gravitational-wave Observations (IndIGO) [63].The increasing number of detectors will boost the capabilities of GW astronomy both in the number of confirmed events (a larger volume of the Universe is cov-ered) and sky localisation (a better triangulation of the source), which will also aid in the search for a counterpart.However, 2G detectors are limited in their sensitivity and future third-generation (3G) ground-based detectors are designed to become more sensitive, precise and capable of probing a larger range of frequencies.Special emphasis should be given to the Einstein Telescope (ET), which is expected to improve the current sensitivity by a factor of 10 [64].ET will also extend the redshift range, e.g.z ∼ 5 for binary black-holes (BBHs) compared to z ∼ 0.5 for 2G detectors [65].The number of detectable multimessenger events is expected to reach tens of thousands of standard sirens [66].While these groundbased detectors will cover a frequency band in the range 1 ≲ f ≲ 10 3 Hz [67], the upcoming space-based 3G detectors, such as the Laser Interferometer Space Antenna (LISA) [68] will have a peak sensitivity near 10 −3 Hz and will be able to detect GW events beyond z = 20, probing a wide range of targets.There are many proposals of 3G GW observatories, such as the DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) [69].However, we have opted to focus our analysis on ET and LISA covering groundand space-based experiments.
In this paper, we aim at forecasting the constraining power of future 3G detectors; in particular we will focus on ET and LISA.Given the potential of such missions we are interested in assessing their ability to constrain modifications to General Relativity (GR) as well as to provide complementary constraints on H 0 using standard sirens.We investigate four models characterised by coupling functions between the DM and DE fields: a conformally coupled quintessence field [21,22], characterised by a conformal coupling in the form of an exponential function of the scalar field; a kinetic model [70] with a conformal function given by a power law of the kinetic term of the scalar field; a purely disformally coupled quintessence field [29,38] with a constant disformal coupling and a mixed disformally coupled quintessence [42,43] which combines the previous model with an exponential conformal coupling.All the scenarios considered are characterised by the same simple exponential potential which introduces one more free parameter.The models we consider differ considerably in the way the (effective) coupling between DM and DE evolves and, in particular, they differ in their background evolution.As such, they represent a well studied sample of models of interacting dark energy suitable for our analysis.We construct our pipeline following the methodology presented in [71][72][73][74] providing for the first time GWs forecasts on the free parameters of the four models in question.
This paper is organised as follows.We start by giving a brief introduction to the physics of standard sirens in Sec.II.Sec.III provides an overview of the methodology used and the details on the simulation of the standard siren events developed for this study, as well as a brief account of the data set combinations considered.We outline the criteria for particular catalogue choices, and discuss the sampling method employed for the forecasts.In Sec.IV we introduce each of the four models under study and present the results of our analysis, emphasising their significant implications.Lastly, in Sec.V, we summarise our results and outline our concluding thoughts and future prospects.

II. GRAVITATIONAL WAVES AS STANDARD SIRENS
Interferometers are sensitive to the strain, h(t) from a GW event, which in the transverse-traceless gauge is described as, [71] where θ 0 and ϕ 0 define the initial location of the event relative to the detector in polar coordinates, ψ is the polarisation of the GW event, and t is cosmic time.We adopt a random sampling method in the range [0−2π] for θ 0 and [0−π] for both ϕ 0 and ψ.The factors F ×,+ describe the antenna beam pattern function, The superscript number indicates which interferometer is being considered, e.g LISA only has two separate interferometers and therefore F (3) = 0. Since the detectors are spatially distributed in an equilateral triangle formation, the other two antenna pattern functions relate to F ×,+ as ×,+ (θ, ϕ + 4π 3 , ψ) . ( As LISA is sensitive to lower frequencies, and equivalently larger masses, it can detect GW events of inspiral mergers lasting over several months, during which the interferometer's position will change relative to the event.This change in position is accounted for following the method described in [72].The timescale of the event is described as Here t c is the time of the merger, t indicates the time at which LISA detects the merger, f is the frequency of the GW, and M c is the chirp mass.The location angles are updated accordingly: which, in turn, are used to update the beam pattern functions.Here we have specified the period, T , as the orbit around the Sun.While the individual masses of the objects are not directly discernible, GW detectors are sensitive to the chirp mass, a collective mass quantity related to the frequency evolution of the signal emitted before the merger, during the inspiral phase of the binary [75], defined as where (1 + z) is a conversion redshift factor from the physical to the observational chirp mass.The Fourier transform of the strain using the stationary phase approximation [73] reads where Ψ(f ) is the phase of the waveform.Notice that when H is inserted into Eq.( 14), the exponential term disappears, meaning that the Ψ(f ) factor can be discarded for this analysis.A is the Fourier amplitude of the waveform, where d L is the luminosity distance from the merger and l is the inclination angle, which we sample randomly between [0 • , 20 • ], as that is the maximum detection inclination range.LISA has been designed to effectively measure frequencies as low as f min = 1×10 −4 Hz, which is why it stands as a promising probe of extreme mass ratio inspiral (EMRI) and binary massive black hole (BMBH) mergers.For the purpose of the simulations, the upper bound frequency of LISA is determined by two quantities: the structure of LISA itself and the last stable orbit of the merging system.LISA can detect frequencies up to f max = c (2πL) −1 , where L is the length of LISA's interferometer arm, taken to be 2.5 Gm and c is the speed of light.Moreover, the total mass of an orbiting system is inversely proportional to its measured frequency, implying that even though massive mergers give rise to large detection amplitudes, the frequency will fall below f min .Therefore if the last stable orbit frequency, f LSO = (6 3/2 2πM obs ) −1 , with M obs being the observed total mass, is found to be lower than f min , we disregard that simulated event.If otherwise it lies between f min and f max then f LSO becomes the new maximum frequency for that event.

III. METHODOLOGY AND DATA SETS
Given the main objective of this study, we create simulated data that forecasts the potential future observations of standard siren events.Specifically, we focus on those that could be detected by ET and LISA.Below, we provide a concise overview of the samples we have generated along with the methodology and the data combinations used.

A. Simulated Cosmology
To simulate GW catalogues from future probes of black hole mergers, the following cosmological quantities are required: the redshift of the merger, z, the value of the Hubble rate at merger, H(z), its comoving and luminosity distance, d c (z) and d L (z) respectively, and the cosmic time between the merger and measurement, t.For this purpose, we resort to the public Einstein-Boltzmann code CLASS code1 [76][77][78], which we extend to accommodate general models of interacting dark energy.This new patch is then used to provide a mock Universe adopting a flat ΛCDM as the fiducial model to simulate the GW data, according to the best-fit cosmological parameters of the Planck 2018 data release [5].These are: the Hubble parameter at present time, H 0 = 67.32km s −1 Mpc −1 , the density of baryons, Ω b h 2 = 0.022383 (with h = H 0 /100) and the density of cold dark matter, Ω c h 2 = 0.12011.Furthermore, we are also interested in the derived quantity Ω 0 m = Ω b + Ω c , which for the fiducial Planck case is Ω 0 m = 0.3144.Provided with the background cosmology, we simulate the merger events to determine the redshiftluminosity relation.First, we generate a redshift distribution of events weighted by a probability distribution.The characteristics of these events, such as the chirp mass, are simulated using a uniform distribution.Although each instance of running the script will yield a different set of simulated data, the resulting conclusions will be unaltered as the fiducial parameters constrain the mock data.Once the mergers have been simulated, we emulate the measurement process from the inspiral, yielding the errors associated with each event.As such, simulated data points are removed if they produce a signal-to-noise ratio below the threshold.

B. Distribution of Simulated Merger Events
ET is designed to probe a range of frequencies, f , similar to that of LIGO, thereby probing merger events of nearby compact objects such as binary neutron stars (BNS) in the mass range of [1,2], [1,2] M ⊙ , and black hole-neutron star binaries (BHNS) in the mass range [3,10], [1,2] M ⊙ , with the [•, •] notation indicating the uniformly distributed mass ranges considered.Advanced LIGO claims a ratio of BHNS to BNS merger events of ∼ 0.03 [79].The redshift probability distribution of these events is proportional to where the comoving distance and the Hubble parameter are taken at various redshifts determined by CLASS.R(z) stands for the merger rate, which, at a linear ap- On the other hand, LISA will target lower frequencies when compared with other proposed 3G detectors, implying sensitivity to events from larger mass binary systems since f ∝ M −1 .Therefore we focus on simulating the detection of events from EMRIs and BMBHs in the ranges [1 − 30], [10 4 − 10 8 ] M ⊙ [80] and [10 4 −10 8 ], [10 4 −10 8 ] M ⊙ [81], respectively.The number of detected BMBH to EMRI events is estimated to follow a 2 : 1 ratio according to the mission's proposal [82,83].
Although in principle LISA will also be able to probe mergers of binary intermediate-mass black holes (IMBHs) and binary compact objects, we opt to discard these from the simulations.This is due to the fact that there is no definitive observational proof of IMBH, and expected events from binary compact objects will only be observed at redshifts z ≈ 3 [84].These events are insignificant since we are interested in the higher range of redshifts for our cosmology.
Considering events involving BMBHs only, the redshift probability distributions are based on the histogram for the L6A2M5N2 mission specification [81] which considers three formation processes of BMBHs.We consider the light seed model (pop III) which attributes the formation of BMBHs to the remnants of population III stars around z = 15 − 20.In [81], two additional scenarios for massive black hole formation were investigated, namely delay and no delay scenarios.These cases involve the collapse of gas in a galactic centre at z = 15 − 20, leading to the formation of a black hole through a heavy seed mechanism with and without a delay between galaxy merger and the merger of the central massive black hole.Further information on these scenarios can be found in Ref. [85].
In our investigation we provide mock data and obtain forecasts for both the delay and no delay cases.However, the analysis reveals that the predicted constraining power from these models shows no actual improvement compared to the pop III case.Consequently, in this paper, we focus solely on the pop III model, as it proves sufficient to forecast the constraining power of LISA.

C. Simulation of Measurements and Errors
To simulate the errors associated with the standard siren catalogue, we follow the methodology of [71][72][73][74].An apparent detection of a GW event is assessed by evaluating the signal-to-noise ratio (SNR), ρ, and only confirmed if ρ > 8.The SNR is defined as where the number labels indicate the interferometer being considered.H has been defined in Eq. ( 10) and S h is the noise power spectral density, an SNR weighting function that accounts for the particular properties of the instruments used.For ET, in particular, S h is designed to follow where, x = f /200 Hz −1 , S 0 = 1.449 × 10 −52 Hz, p 1 = −4.05,p 2 = −0.69, a 1 = 185.62,a 2 = 232.56,b n = {31.18,−64.72, 52.24, −42.16, 10.17, 11.53}, and c m = {13.58,−36.46, 18.56, 27.43}, assuming a lower cutoff at f = 1 Hz.On the other hand, for LISA, S h depends on the instrumental (or short) noise, S inst , the noise from low-level acceleration, S acc , and the confusion background noise, S conf [85]: where S acc = 9 × 10 −30 /(2πf ) 4 (1 + 10 −4 /f ), S inst = 2.22 × 10 −23 and S conf = 2.65 × 10 −23 .Therefore, the total SNR contribution for each detector is given by combining (10) with either Eq. ( 15) or Eq. ( 16) for the ET and LISA, respectively: The instrumental error in the luminosity distance is determined via the Fisher matrix, following [72].Since H ∝ d −1 L this results simply in where the factor of 2 accounts for the symmetry in the inclination angle, which actually ranges from −20 • to 20 • .The error due to gravitational lensing is, reduced by a half to account for both the merger and ringdown of the event.
Being space-based, LISA is also subject to an error associated with the peculiar velocities of GW sources [86]: with an estimate of the peculiar velocity of the host galaxy with respect to the Hubble flow of ⟨v2 ⟩ = 500 km s −1 .Bringing all the contributions together, the total error in the luminosity distance is simply a combination of the errors in Eqs. ( 19)-( 21): The simulation allows us to interpolate any number of events over a continuous redshift distribution in the range 0 < z ≲ 5 for ET and 0 < z ≲ 10 for LISA.However, the number of mergers detected by ET will depend on factors such as running costs and the complementary detection with other experiments [71].ET is expected to report more than 10 4 mergers yearly.However, due to the scarcity of EM counterpart signals, the predicted number of detectable mergers with an actual EM counterpart over the course of 10 years is approximately 200 [87].According to [81], LISA's number of detected mergers, for a 10-year mission proposal, is 56 events.
To incorporate uncertainty into the luminosity distance of each merger, we apply a Gaussian distribution centred around the background cosmology.The standard deviation for this distribution is set to the calculated errors, σ d L .This introduces artificial randomness around each merger, leading to a larger deviation from ΛCDM in LISA compared to ET.The reason for this difference is that LISA probes larger redshifts, which are associated with larger errors, resulting in a broader spread of the data, as depicted in Fig. 1.

D. Data Sets and Likelihoods
To examine the fit of the simulated data to the coupled quintessence models considered in this study, we employ the Markov Chain Monte Carlo (MCMC) method, using samples generated from our modified version of CLASS interfaced with the Monte Python sampler [88,89].In particular, we resort to the Nested Sampling algorithm through the MultiNest 2 [90-92] and PyMultiNest3 [93] packages to estimate observational constraints on the free parameters, instead of the traditional Metropolis-Hastings algorithm.The Metropolis-Hastings algorithm struggles to explore the full line of degeneracies between the parameters, resulting in false peaks in the posterior distribution, which the sampler cannot move away from.Nested sampling is able to explore the full extent of the degeneracies as it is much better suited for multi-model sampling (see section IV) and other more complicated distributions.Subsequently, we analyse the MCMC chains and present the results using the GetDist4 Python package [94].
The likelihood function for the simulated data set of standard siren GW events is constructed according to the effective Gaussian distribution: where SS (z) is the observed luminosity distance, which in this case corresponds to the samples generated according to the procedure outlined above; d SS (z) is the model-dependent theoretical prediction for the luminosity distance of the event, computed numerically with the modified CLASS code; σ d L is the total error in the luminosity distance, as defined in Eq. ( 22); and n is the number of observed events.
Since we want to forecast the constraining power of standard siren data probed by ET and LISA on coupled quintessence models, we assess the independent and combined constraints with current background data.This allows for a direct comparison of whether GW catalogues will improve the constraints on {Ω 0 m , H 0 } and on the model-specific parameters, affecting the background evolution.In particular, we include baryonic acoustic oscillations (BAO) data from the Sloan Digital Sky Survey (SDSS) DR7 Main Galaxy Sample [95], SDSS DR12 consensus release [96] and the 6dF Galaxy Survey [97], in combination with distance moduli measurements of 1048 type Ia Supernova (SNIa) data from Pantheon [98].This combined data set is referred to as "SNIa+BAO".
Our analysis involves a set of free sampling parameters, including the baseline ΛCDM cosmological parameters (Ω 0 m , H 0 ) and the parameters associated with each coupled quintessence model, for which we consider as fiducial value their ΛCDM limit.The models discussed in Sec.IV reduce to ΛCDM in the following limits: λ = 0 and β = 0 for IV A; λ = 0 and α = 0 for IV B; λ = 0 and D 0 = 0 for IV C; λ = 0, β = 0 and D 0 = 0 for IV D. We adopt flat priors for all parameters within the ranges specified in Table I.

IV. FORECAST RESULTS
In what follows, we employ the methodology and the data sets discussed in Sec.III to investigate the power that LISA and ET standard sirens have in constraining the cosmological parameters {Ω 0 m , H 0 }, the In each subsection that follows, we show the resulting 2D contours, and 1D marginalised posterior distributions for {H 0 , Ω 0 m } plus the set of model-specific parameters (see Table I) for the cases of ET and LISA and their combination.These plots also include a combination of SNIa+BAO, and for reference, the results of SNIa+BAO alone.The results are also summarised in a table with their corresponding 1σ, identified in the text with the notation {σ p } where p is an index spanning over the model parameters.We also use p } where here i and j stand for two different data sets, to denote the change in error for the specific parameter p.

A. Conformal Coupling
The first model we consider is the conformal coupling model, for which where C(ϕ) is defined according to Eq. ( 1) and V (ϕ) is the DE potential energy.The exponential parameters β and λ are constant dimensionless parameters and V 0 is a constant with dimensions of (mass) 4 that sets the energy scale of the potential 5 .In such models, the mass of the DM particles becomes ϕ-dependent and the DE field mediates a long-range force between DM particles so that the effective gravitational coupling is given by 99,100].The free parameters we are particularly interested in are the slope of the potential λ and the coupling parameter β.Constraints on this model have been obtained in Ref. [101] using background data only (H(z), BAO and supernova Union2.1).Using these data, the authors found the following upper limits: β < 0.193 and λ < 1.27.In [102] stronger constraints have been obtained using Planck data, BAO and SNIa data, also in line with Ref. [103], in which the authors found β < 0.0298 and λ < 0.6 for the 1σ upper limits.
Compared to ET+SNIa+BAO and LISA+SNIa+BAO to SNIa+BAO, there is also a reduction in σ for all parameters except one.The ET+SNIa+BAO data set results in a nominal increase in σ Ω 0 m compared to SNIa+BAO.Comparing the errors of ET and LISA to SNIa+BAO we see only minor changes in the constraining power regarding Ω 0 m , with ET performing slightly worse and LISA slightly better.A similar trend occurs for the parameter λ, with ET performing nominally worse and LISA better.However, particular attention should be given to the significant reduction in σ H0 when comparing ET and LISA to SNIa+BAO.There is a reduc-tion in the error by a factor of F (SNIa+BAO,ET) H0 = 0.090 and F (SNIa+BAO,LISA) H0 = 0.11.Forecasting GWs will improve the constraints on H 0 , suggesting that GWs will be critical in addressing the Hubble tension.On the other hand, we see the opposite effect with β, with an increase in the error by a factor of F = 0.32.In comparing the constraining power of ET and LISA (see Fig. 4), it is evident that they have comparable spreads for the cosmological parameters.An interesting feature we observe is that ET is more constraining in regard to H 0 .We attribute this feature to the fact that the ET catalogue has more data points than LISA at low redshifts, as illustrated in Fig. 1.
By combining GW data from LISA and ET, which implies an increase of data points over a wide range of redshifts, we predict an enhanced constraining power in the cosmological parameters, {H 0 , Ω 0 m }, compared to the SNIa+BAO case, more precisely F (SNIa+BAO,ET+LISA) Ω 0 m ,H0 = {0.65,0.063}.However, for the model parameters β and λ, we observe modifications to the constraining power with F (SNIa+BAO,ET+LISA) β,λ = {1.8,0.86}.The combination of ET+LISA with SNIa+BAO results in a negligible change of constraining power for Ω 0 m , H 0 , λ.Only the parameter β is more constrained when the data sets are combined, with the 1σ reduced by almost a third.

B. Kinetic Conformal Coupling
As an example of a coupled quintessence model in which the conformal function is less trivial, we focus on a pure dependence on derivatives of the scalar field through the kinetic term of ϕ, X = −∂ µ ϕ∂ µ ϕ/2, to which we refer as the kinetic coupling.Such a setting has been proposed in [104] (see references therein as well), and we focus on the particular example of a power law, as studied in [70].Even though this model is proposed based on a Lagrangian framework (L DM −→ X/M 4 Pl α L DM ), at the background level it is equivalent to the kinetic-dependent conformal transformation ḡµν = C(X)g µν , with where α is a dimensionless constant and a simple exponential potential has been assumed just like in the previous case, and the same considerations apply for λ and V 0 .
In summary, an analysis based on Planck and the SNIa+BAO background data in Ref. [70] reveals the power of BAO data in constraining Ω 0 m , which is highly correlated with the steepness of the potential λ.The coupling parameter α is constrained to be of the order of 10 −4 .The constraints on the cosmological parameters are found to be compatible with the ΛCDM ones within the errors.Moreover, a positive correlation between H 0 and Ω 0 m is identified.While this trend is attributed to the evolution of the linear perturbations for non-vanishing α, we find that it is still present for the background standard siren data sets.
From the results presented in Figures 5, 6 and 7, and summarised in Table III, we analyse the constraints on the parameters {Ω 0 m , H 0 , λ, 10 4 α} for the same data sets as in the previous case.When evaluating the errors from ET standard sirens and comparing them to SNIa+BAO data, we observe that for most parameters, ET's 1σ constraints are of the same order, apart from the H 0 parameter, which is improved by 1 order of magnitude.This reduction is quantified by the fractional change of F (SNIa+BAO,ET) Ω 0 m ,H0,β,λ = {1.1,0.12, 1.0, 1.1}.When the data sets are combined (ET+SNIa+BAO), we find that the 1σ region is narrower for all parameters compared to ET alone, with F (ET, ET+SNIa+BAO) Ω 0 m ,H0,α,λ = {0.84,0.92, 1.0, 0.92}.In the case of LISA standard sirens, we observe that all cosmological and model parameters are better or equally constrained by LISA alone compared to SNIa+BAO, with F (SNIa+BAO,LISA) Ω 0 m ,H0,α,λ = {0.91,0.13, 1.0, 1.0}.Combining LISA with SNIa+BAO, we find improved constraints with respect to the SNIa+BAO data set alone.Moreover, when comparing LISA+SNIa+BAO with LISA alone, the former shows an even better constraining power, with the most significant reduction in error observed for Ω 0 m , with F (LISA,LISA+SNIa+BAO) Ω 0 m ,H0,α,λ = {0.78,0.92, 1.0, 0.87}.For both the ET and LISA data sets, the accuracy on H 0 can be improved by 1 order of magnitude (0.36 for ET and 0.39 for LISA) compared to SNIa+BAO (3.1), as reported in section IV A as well.Interestingly, for all data sets and combinations, the accuracy of the model parameters remains largely unaffected, with the 1σ region for λ showing only nominal changes and remaining unchanged for α.Comparing the constraining power of ET and LISA with their combination, ET+LISA, we see that the latter provides better constraining power for the cosmological parameters than any of the other data sets analysed.Regarding the model parameters, there seems to be a minimal change in accuracy compared to the single ET or LISA data sets.We do note that the GW combination provides better accuracy with respect to SNIa+BAO, with F (SNIa+BAO,ET+LISA) Ω 0 m ,H0,α,λ = {0.68,0.084, 1.0, 0.91}.The full combination of ET+LISA+SNIa+BAO, has a negligible change in the constraints when compared to ET+LISA for all parameters.
Comparing the accuracy of the constraints of the Kinetic model obtained in Ref. [70] with CMB TT, TE and EE Planck 2018, Planck CMB lensing, BAO and SNIa data, we note that the parameter α is better constrained by CMB data and its combination with BAO and SNIa by 1 order of magnitude when compared to our data combinations, given that the latter only depend on the background evolution.More precisely, we report σ 10 4 α = 2.9 for all the data set combinations while in Ref. [70] this was reduced to σ 10 4 α = 0.95 (Plk18), 0.84 (Plk18+SNIa+BAO), and 0.7 (Plk18+SNIa+BAO+Lensing).Future ET and LISA catalogues will be able to constrain λ at the same level as Planck CMB data (σ λ = 0.48 with Plk18 and σ λ = 0.2 with both Plk18+SNIa+BAO and Plk18+SNIa+BAO+Lensing). On the other hand, the standard siren data will better constrain H 0 by 1 order of magnitude with respect to Plk18 (σ H0 = 2.5).CMB lensing data increase the constraint by 1 order of magnitude, namely with accuracy σ H0 = 0.6, which is of the same order of magnitude as the ET and LISA cases, even though the standard sirens perform better in terms of the relative error with σ H0 < 0.4 for all the combinations considered.

C. Disformal Coupling
In the following we study the model with disformal coupling only, (26) in which case the conformal contribution vanishes and D is simply a constant with dimensions of (mass) −4 in Eq. ( 2) (and hence D 0 has units of (mass) −1 ), and V (ϕ) follows the same considerations as in the previous cases.The constraints on this model have been obtained in [101] and [102].It was found that using background data only (H(z), BAO and supernova Union2.1 data) results in the following constraints: D 0 > 0.07 meV −1 and λ < 1.56 at 95.4% [101].An upper bound can be obtained for D 0 with CMB data (including lensing) and BAO, SNIa, cosmic chronometers, cluster abundance, and H 0 priors which is D 0 < 0.2500 meV −1 and a stringent upper limit for λ is < 0.6720 at 1σ [102].
From Figs. 8, 9 and 10, summarised in Table IV, we analyse the results for the parameters {Ω 0 m , H 0 , D 0 , λ} for the same data sets as in the previous cases.In our analysis of the ET data set alone, we observe an improved accuracy for all cosmological and model parameters, {Ω 0 m , H 0 , D 0 , λ}, compared to SNIa+BAO, with F (SNIa+BAO,ET) Ω 0 m ,H0,D0,λ = {0.71,0.10, 0.98, 0.85}.As expected, the combination of data sets, ET+SNIa+BAO, also results in improved accuracy compared to SNIa+BAO.Compared to the ET data set alone, there are only minor changes in the parameters' accuracy, F In the case of LISA standard sirens, we find that the cosmological and model parameters follow a similar accuracy trend with F (SNIa+BAO,LISA) Ω 0 m ,H0,D0,λ = {0.71,0.11, 0.98, 0.85}.Moreover, the same is true for the combination LISA+SNIa+BAO, with F (LISA,LISA+SNIa+BAO) Ω 0 m ,H0,D0,λ = {1.0,1.0, 0.98, 1.22}.There is only a nominal change in the accuracy of parameters compared to that of LISA alone, apart from λ, which results in a larger 1σ region, with σ λ = 0.71.Regardless of the data combination, {Ω 0 m , D 0 , λ} are constrained at the same level, with the parameter λ having a slight improvement in accuracy for both ET and LISA (both have σ λ = 0.58 compared to σ λ = 0.7 for SNIa+BAO).The accuracy of the H 0 parameter is 1 order of magnitude better for both ET and LISA than SNIa+BAO.
There is no change in the model parameters for ET and LISA, and thus we see no noticeable change in the constraints for ET+LISA.However, there is an increase in accuracy for the cosmological parameters.As both ET and LISA improved the constraints compared to SNIa+BAO, we observe the expected result, that ET+LISA have further improved constraints with F (SNIa+BAO,ET+LISA) Ω 0 m ,H0,D0,λ = {0.55,0.071, 0.96, 0.85}.Following the trend with LISA and the combination with SNIa+BAO, we note that the combined data set ET+LISA+SNIa+BAO, has very little change in the accuracy compared to ET+LISA apart from the constraint for λ, which results in a worse accuracy than ET+LISA and SNIa+BAO with F (ET+LISA,ET+LISA+SNIa+BAO) Ω 0 m ,H0,D0,λ = {1.0,1.1, 1.0, 1.3}.
The first thing to be noted in comparison with the results reported in Ref. [101] is that we were able to derive constraints at 68% C.L. and not only at 95% C.L. for all the model parameters, therefore providing better constraints in all the cases.Moreover, the constraining power on H 0 is largely improved from σ H0 ≈ 2.2 for the background data to σ H0 ≈ 0.3 in all the cases, including standard sirens.When compared with the results of CMB, CMB lensing and additional data in Ref. [102], which report only upper bounds for λ and D 0 , we see that both parameters are constrained at 68% C.L. with standard sirens, with lower and upper bounds, in particular with more accommodating upper bounds, as this analysis includes only back-  ground data.Accordingly, the error in H 0 is brought to the same order of magnitude with σ H0 ≈ 0.9, which is still about 3 times larger than the one reported in Finally, we discuss a model with a mixed coupling consisting of a conformal and a disformal part.Specifically, we consider As for the disformal case, the constraints on such a model were discussed in [101] and [102].For the same background data it has been reported that D 0 > 0.102 meV −1 , β < 0.453 and λ < 1.59 at 95.4% [101].An example of the constraints including CMB data in [102] are β ≲ 0.17 and λ ≲ 0.35 at 1σ, with the details depending on the data sets used, with the disformal coupling D 0 not always well constrained for this case, with lower bounds of D 0 ≳ 0.35 meV −1 for some data combinations.
In Figs.11, 12 and 13 and Table V we show the results for the parameters {Ω 0 m , H 0 , β, D 0 , λ} for the same data sets as before.In our analysis of the ET data set alone, we observe improved accuracy for the cosmological parameters, Ω 0 m and H 0 , compared to the SNIa+BAO data set, with F (SNIa+BAO,ET) Ω 0 m ,H0 = {0.61,0.088}.The combined data sets, ET+SNIa+BAO, show comparable results, with a slight increase in accuracy compared to ET alone.For the model parameters, {β, D 0 , λ}, ET compared with SNIa+BAO demonstrates close constraining power, with F (SNIa+BAO,ET) β,D0,λ = {1.0,0.92, 1.0}.However, the combined data set leads to an increase in the error in β with σ β = 0.71 for ET+SNIa+BAO.For the case of LISA standard sirens, we find that the cosmological parameters follow a similar trend as ET, with increased accuracy, F (SNIa+BAO,LISA) Ω 0 m ,H0 = {0.72,0.088}, with the combined data set showing a comparable trend.However, it is worth noting a noticeable reduction in accuracy for H 0 between LISA and the combined data set, F (LISA, LISA+SNIa+BAO) H0 = {1.1}.Regarding the model parameters, {β, D 0 , λ}, we find that unlike for ET, LISA alone exhibits increased accuracy compared to SNIa+BAO, except for D 0 , F Regarding the model parameters, we find only a very small change compared to SNIa+BAO, F (SNIa+BAO, ET+LISA) β,D0,λ = {0.96,1.2, 1.1}, noting that both D 0 and λ are slightly less constrained.
In summary, regardless of the combination of data sets considered, the constraints on Ω 0 m , β, D 0 , λ are of the same order of magnitude as those obtained from SNIa+BAO.Additionally, we note that the accuracy on H 0 is improved by 1 order of magnitude for both ET and LISA compared to SNIa+BAO.
Similarly to the comparison in Sec.IV C, the main improvement in contrast with the results reported in Ref. [101] is the fact that we can obtain constraints at 68% C.L. for all the model parameters.The potential parameter λ is constrained with upper bounds for the standard sirens at 1σ of the same order of the 2σ ones reported in the previous studies.Moreover, the constraining power on H 0 is largely improved from σ H0 ≈ 2.1 for the background data to σ H0 ≈ 0.3 in all the cases, including standard sirens.The comparison with results including CMB, CMB lensing and additional data in Ref. [102], which are either unable to constrain D 0 or find just a lower bound and report only upper bounds for λ and β, shows that standard sirens successfully constrain the three model parameters at 68% C.L. for all the combinations, which is a great improvement given that only background data has been considered.Including CMB data brings the error in H 0 to the same order of magnitude with σ H0 ≈ 0.6, which is still around 2 times larger than the ones reported in this analysis.

V. CONCLUSIONS
In this paper, we have explored the potential of future GWs detectors, namely LISA and ET, to constrain conformal and disformal couplings between the dark energy and dark matter fluids.We have considered four models: conformal coupled quintessence, a kinetic model, constant disformal coupled quintessence and a mixed conformal-disformal model.All the cases considered have the same exponential potential of which we have constrained its slope, λ.
We have generated mock catalogues of standard siren events with ET and LISA specifications and with those we performed an MCMC analysis con- sidering, separately, their combination with current SNIa+BAO data as well for reference.Under the assumption we have used to generate the mock data assuming a particular cosmology, we find the following: • The conformal coupled quintessence: the combinations of LISA+SNIa+BAO and ET+SNIa+BAO improve the constraints on both the slope parameter, λ, and the conformal coupling parameter β.The combination of ET+LISA with SNIa+BAO reduces the error in β by one third.
• The kinetic model: ET and LISA alone cannot improve the constraints on λ or on the conformal all combinations can constrain the disformal parameter D 0 at 1σ with the same order of magnitude and a small improvement for LISA+SNIa+BAO.For the slope parameter instead, both ET, LISA and their combination perform better than SNIa+BAO.Moreover for the full catalogue combination of ET+LISA, the error in Ω 0 m can be reduced.• The mixed conformal-disformal coupled quintessence: for all the parameters of the model there is not a significant improvement on the accuracy of their constraints when using ET or LISA data separately.A small reduction in the size of the 1σ region is reported only for the disformal parameter, D 0 , in the full combinations.The error in Ω 0 m is slightly reduced for ET+LISA.
Regardless of the model considered we found that the accuracy on the H 0 parameter increases by 1 order of magnitude at 1σ when compared to the combination of BAO and SNIa data.This is promising in light of solving/understanding the H 0 tension.This improvement is also responsible for the increased accuracy on the constraints for the model parameters when we consider the full combinations that we just reviewed.
Ultimately, our results show that future 3G detectors can improve our knowledge on DE-DM interaction and shed light on the H 0 tension.
0. Nonetheless, when the background data is combined with ET and/or LISA, the constraints improve by F

FIG. 3 .
FIG. 3.68% and 95% C.L. 2D contours and 1D marginalised posterior distributions for the parameters {H0, Ω 0 m , λ, β} in the conformal coupled quintessence model with LISA mock data (charcoal filled line), SNIa+BAO data (red dotted line) and their combination (yellow dashed line).The scale is the same as in Fig.2for comparison purposes, with the SNIa+BAO contours standing as the reference.The dotted lines depict the fiducial values for the mock data {Ω 0 m , H0} = {0.3144,67.32}.

(
SNIa+BAO,LISA) β,D0,λ = {0.98,1.1, 0.98}.The combination of the data sets results in comparable accuracy to LISA alone, with F (LISA,LISA+SNIa+BAO) β,D0,λ = {1.0,0.89, 1.0}.Similar to sections IV A to IV C, the combination of the GW data sets leads to a significant change in the accuracy of Ω 0 m and H 0 compared to the SNIa+BAO data sets, F (SNIa+BAO, ET+LISA) Ω 0 m ,H0 = {0.49,0.061}.The accuracy of the cosmological parameters is also slightly enhanced compared to ET or LISA alone.

TABLE I .
Flat priors on the cosmological and model parameters sampled in Sec.IV.model-dependent conformal and disformal coupling parameters and the steepness of the self-interacting potential.In particular, we consider four interacting DE models: a standard coupled quintessence model, a kinetically coupled model, a constant disformal model and a mixed conformal-disformal model.For each of the four scenarios, we provide a brief review of the theoretical framework before presenting the forecasts obtained considering the specifications and assumptions discussed in previous sections.

TABLE II .
Marginalised constraints on cosmological and model parameters for the Conformal Coupled Quintessence model at 68% C.L.

TABLE III .
Marginalised constraints on cosmological and model parameters for the Kinetic Model at 68% C.L.

TABLE IV .
Marginalised constraints on cosmological and model parameters for the Constant Disformal Coupled Quintessence Model at 68% C.L.

TABLE V .
Marginalised constraints on cosmological and model parameters for the Mixed Conformal-Disformal Coupled Quintessence Model at 68% C.L.