Cosmological perturbations in the energy-momentum squared gravity theory: constraints from gravitational wave standard sirens and redshift space distortions

Recently, the energy-momentum squared gravity (EMSG) has emerged as a novel approach for addressing the problem of the Big-Bang singularity by introducing the self-coupling of matter. This theory introduces a scalar, denoted as $ \mathbb{T}\equiv T_{\mu\nu}T^{\mu\nu} $, which is constructed from the energy-momentum tensor (EMT) and distinct from theories that rely on the trace of $ T_{\mu\nu} $. Its modifications to general relativity (GR) result from the matter instead of the geometry, which results in interesting cosmological outputs. Reference [1] demonstrated that a simple model of the EMSG theory in the metric formulation yields a viable cosmic behavior, presenting a genuine sequence of cosmic epochs. Notably, the key deviation from the standard cosmological picture lies in the absence of an early universe singularity. In essence, the EMSG theory posits the existence of a minimum length and a maximum, yet finite, energy density during the early universe. Presently, interest in this theory is not limited to solving the Big-Bang singularity. Depending on specific EMSG model under consideration, it may propose intriguing modifications to the entire cosmic history, extending beyond the early universe. While the EMSG theory does not evade black hole singularities, it predicts larger masses for neutron stars under ordinary equations of state [2, 3]. For recent work on quark compact stars within the framework of the EMSG theory, readers are referred to Ref. [4].

Various aspects of the EMSG theory have been explored. In Ref. [5], different EMSG models were scrutinized using the dynamical system approach, revealing intriguing cosmological behaviors through the analysis of relevant fixed points. Additional cosmological investigations related to the EMSG theory are available in Refs. [6−12]. Reference [13] investigated in depth the Jeans analysis within the context of the EMSG theory, introducing a novel Jeans mass. The study also explored the local stability of hyper-massive neutron stars by identifying a generalized version of Toomre's parameter in the EMSG theory, outlining constraints on the theory's free parameters. A notable finding reported in Ref. [12] indicates that the bounce mentioned in Ref. [1] does not accurately describe the actual universe. Specifically, the EMSG theory exhibits a viable bounce only in a matter-dominated cosmological toy model [12]. To achieve a viable bounce, we should replace the quadratic scalar $ \mathbb{T} $ in the action with the term $ \mathbb{T}^{5/8} $. Reference [14] investigated the motion of light in the weak-field limit of the EMSG theory and showed that this theory passes the solar system tests with flying colors. The constraints on the EMSG theory by binary pulsar observations was studied in Ref. [15]. Moreover, Ref. [16] presented an analysis on the gravitational-wave (GW) radiation and radiative behavior of relativistic compact binary systems in the so-called scale-independent EMSG theory. Please see Refs. [17−25] for more investigations on the EMSG theory.

Moreover, the linear perturbation is an effective method for discriminating different modified gravitational theories as they may have the same cosmological background solution, but the evolutions of the perturbed modes are different. Generally, the linear perturbation can be decomposed into scalar, vector, and tensor perturbations. Because the vector mode decays with the cosmological expansion, most studies focus on the scalar and tensor cosmological perturbations. Tensor perturbations correspond to the propagation of GWs in the cosmological background, serving as standard sirens for determining the distance-redshift relation [26]. Therefore, with the upcoming deployment of ground- and space-based GW detectors such as the Einstein Telescope (ET) [27], the Laser Interferometer Space Antenna [28, 29], the DECI-hertz Interferometer Gravitational-wave Observatory [30], the Big-Bang Observer [31], Taiji [32−35], and TianQin [36−42], we anticipate the capability to measure cosmological expansion with high precision up to significant redshifts. Furthermore, the growth rate of matter density perturbation provides an effective approach for estimating the distribution of matter in the universe [43] and is crucial for theoretically differentiating between various gravitational theories [44−47]. The growth of matter density perturbation is highly dependent on the underlying theory of gravity, whether it is GR or a modification of it. Modified field equations result in distinct evolutions of cosmological perturbations, altering the characteristic imprints left in the cosmic microwave background (CMB) and in the matter power spectrum inferred from galaxy clusters. Therefore, it is distinguishable as an exceptional probe for modified gravitational models, facilitating discrimination between these models and conventional dark energy models.

Numerous studies have extensively explored the viable cosmological solutions within EMSG theories, and the constraints on these models with the cosmological electromagnetic observations, such as the cosmic chronometer [8], CMB, baryon acoustic oscillations [9], and type Ia supernovae [10], have been investigated. However, the constraints from GW signals on the EMSG theories have not been investigated yet and it would be interesting to know whether the propagation of GWs is affected. Moreover, the modifications on the matter stress affects the growth of matter perturbations. Thus, constraining the EMSG theories using the redshift space distortion (RSD) observations is effective.

In this paper, we aim to derive the propagation of GWs and the linear evolution for the matter density perturbation and then proceed to extract the observational constraints on the model parameters of some explicit EMSG models, utilizing the mock standard sirens data and the actual RSD data. Moerover, although the spherical overdensities within the EMSG theory was investigated in Ref. [48], the differential equation for the matter density perturbation is derived from the usual continuity equation, which is only valid for some special cases of the EMSG models. In this paper, we begin with the non-conserved EMT and obtain the general equations of motion for the matter density perturbation.

The remainder of this paper is organized as follows. In Sec. II, we briefly review the EMSG theory and present the solutions of two power-law EMSG models. In Sec. III, we derive the tensor and matter density perturbations of the theory. In Sec. IV, we present the formalism to simulate GW detection from the ET and the method for RSD data analysis; subsequently, we determine the fit values of the model parameters by using the mock GW and actual RSD data. Section V presents the conclusions.

Let us start with the action for the EMSG theory [1, 6]:

where $ \mathbb{T}\equiv T_{\mu\nu}T^{\mu\nu} $, and $ T_{\mu\nu} $ is the EMT, which represents the same matter with the Lagrangian $ \mathcal{L}_m $. Varying the action with respect to the metric $ g^{\mu\nu} $, we can derive the field equations as

where $ f_R\equiv\dfrac{\partial f}{\partial R} $, $ f_{\mathbb{T}}\equiv\dfrac{\partial f}{\partial \mathbb{T}} $, and

The matter Lagrangian giving the perfect-fluid EMT is not unique as both of $ \mathcal{L}_m=p $ and $ \mathcal{L}_m=-\rho $ can provide the same EMT. In this paper, we consider $ \mathcal{L}_m = p $ for consistency. Thus, the tensor $ \theta_{\mu\nu} $ reduces to

The covariant divergence of Eq. (2) is

The nonvanishing of the right-hand side of the above equation indicates explicitly that the covariant conservation of the EMT is not a priori in the EMSG theory. These nonvanishing terms act as source terms accounting for the stress-energy transfer between matter and geometry. To proceed, we must specify explicit models of this theory. Without loss of generality, we will focus on the power-law one without the mixed terms between R and $ \mathbb{T} $, which is one of the most viable models for cosmology [1, 6−10, 12, 48], although our methodology can be applied for any functional form of $ f(R, \mathbb{T}) $. Additionally, the matter density perturbation will be considerably simplified without the mixed terms. The power-law scenario corresponds to [7, 8]

where $ \kappa=8\pi G $ with G being Newton's constant, η is the coupling constant signifying the strength of the EMT-powered modification to gravity, and n is the power of $ \mathbb{T} $. Generally, this system can only be solved numerically. However, analytical solutions can be determined for some particular values of n, i.e., $ n=0 $ [8], $ n=1/2 $ [6], and $ n=1 $ [1]. Clearly, for $ n=0 $, this theory reduces to GR, which is the trivial case. The scenario with $ n=1/2 $ corresponds to the so-called scale-independent EMSG theory [9], which exhibits the potential to be influential across various cosmological epochs characterized by different energy density scales and plays a significant role in both early and late-time evolutions of the universe. The model with $ n=1 $ corresponds to the original EMSG theory, which can pass the solar-system weak-field tests [14, 49]. In the following, we will take $ n=1/2 $ and $ n=1 $ as examples to investigate the linear cosmological perturbations of the EMSG theory.

For $ n=1 $, the covariant conservation of the EMT can be accomplished, and the usual continuity equation is satisfied, which leads to the solution for the matter density $ \rho=\rho_0 a^{-3} $. The Friedmann equation in terms of the redshift can be expressed as [1]

where $ \Omega_m\equiv\dfrac{\kappa \rho_0}{3H_0^2} $, $ \Omega_{\Lambda}\equiv\dfrac{\Lambda}{3H_0^2} $, and $ \Omega_Q\equiv\dfrac{\eta \rho_0^2}{6H_0^2}=\dfrac{3\eta H_0^2\Omega_m^2}{2\kappa^2} $, which satisfy $ \Omega_m+\Omega_Q+\Omega_{\Lambda}=1 $.

For $ n=1/2 $, Eq. (5) leads to $ \rho=\rho_0 a^{-3(1+\eta_t)} $ with $ \eta_t\equiv\dfrac{\eta}{2\kappa} $, and the Friedmann equation is given by [6]

In this section, we establish the key equations governing the linear cosmological perturbations of the gravitational field in the EMSG theory. These linear perturbations will be assumed to evolve on a homogeneous and isotropic background described by the Friedmann-Lemaître-Robertson-Walker metric. According to the $ S O(3) $ symmetry of the background, these small perturbations can be decomposed into scalar, vector, and tensor perturbations, which can be studied independently. Because the vector mode decays with the cosmological expansion, we primarily restrict our attention to the tensor and scalar modes of the cosmological perturbations.

Let us focus on the equations of motion of GW first. The perturbed metric can be written as

where $ h_{ij} $ represents the GW perturbation, which satisfies the transverse-traceless conditions, i.e., $ \eta^{ij}h_{ij}=0 $ and $ \partial^i h_{ij}=0 $. Subsequently, the field equations for $ h_{ij} $ derived from Eq. (2) is

where the dimensionless parameter β is defined as $ \beta\equiv -\dfrac{\dot{f}_R}{3H f_R} $, and the dot indicates the derivative with respect to t. Equation (10) can be rewritten in terms of the conformal time as follows:

where τ is the conformal time defined by ${\rm d}\tau={\rm d}t/a(t)$. The above equation indicates that the speed of a GW is equal to the speed of light, and thus the constraints from the observation of GW170817/GRB170817A is trivially satisfied in the EMSG theory [50]. Because the presence of the friction term, the GW luminosity distance is given by [51, 52]

Here, $ d_L^{\rm EM} $ is the standard luminosity distance for an electromagnetic signal, namely

which is the same for the gravitational radiation in GR.

Next, let us consider the evolution of the matter density perturbation for the EMSG theory. In the Newtonian gauge, the metric of the linear scalar perturbations can be expressed as

where $ \Phi(t, x_i) $ and $ \Psi(t, x_i) $ represent the scalar perturbations. The components of the perturbed EMT can be given by

where we have assumed that the universe at the time in which we are performing the perturbations is in the dust dominated phase. By considering Eqs. (2) and (14), we can obtain the perturbed equations with the components $ (00),\; (ij),\; (ii) $, and $ (0i)=(i0) $, where $ i,\; j=1,2,3 $, $ i\neq j $ in Fourier space, respectively given by

with

We emphasize again that the above perturbed equations ignores the mixed terms between R and $ \mathbb{T} $ because we are focusing on the power-law model. Moreover, we can easily observe that by taking $ f(R, \mathbb{T})=f(R) $ in the above results, the perturbed equations in $ f(R) $ gravity are recovered [53]. The covariant divergence of the EMT renders to the following first-order equations:

for the temporal and spatial components, respectively.

From Eqs. (17), (18), and (21), and considering the subhorizon approximation, which is given by $ H^2\Psi,\; H\dot{\Psi},\; \ddot{\Psi}\ll k^2\Psi/a^2 $, we can obtain

Subsequently, combing Eqs. (22) and (23) and inserting the last two expressions, we can derive the equation of motion for the matter density perturbation ($ \delta_m\equiv \delta\rho/\rho $) as

where $ A(t) $ and $ B(t) $ are presented in the Appendix, and this differential equation can be solved by setting the initial conditions in the deep matter dominated era [54, 55]. In the following, we will transform this equation to redshift coordinates and analyze the dynamics of the matter density perturbation as a function of the redshift.

To compare the model with the observational data, we use the data on $ f\sigma_8 $ [56−58], defined as

where $ \sigma_8\equiv\sigma_8(0) $ denotes the current value. The comoving wavenumber is set to $ k=0.1h $ Mpc$ ^{-1} $ to ensure that we are dealing with the linear regime [59, 60], and h represents the reduced Hubble constant satisfying $ H_0=100 h $ km/s/Mpc. In the following, we will estimate the fit values of the parameters $ H_0 $, $ \Omega_{m}\equiv \Omega_{m}(0) $, $ \sigma_8 $, and $ \eta_t $ or $ \Omega_Q $ with the aid of the simulated GW and observed RSD data.

In this section, we first briefly review the simulation of GW detections from the ET. Thereafter, we present the settings for the RSD data analysis. Finally, we explore the constraints on the model parameters through a Markov Chain Monte Carlo (MCMC) analysis using the data from the simulated GW detections and the observed RSD datasets.

Following Refs. [61, 62], the masses of a neutron star (NS) and black hole (BH) are assumed to be uniformly distributed in the ranges $ [1,2] M_{\odot} $ and $ [3,10] M_{\odot} $, respectively. Here, $ M_{\odot} $ is the solar mass. The ratio between NS-NS and BH-NS events is taken to be $ 0.03 $ according to the prediction of the Advanced LIGO-Virgo network [63]. The redshift distribution of the sources assumes the form

where $ d_C $ denotes the comoving distance, which is defined as $d_C(z)\equiv \int_0^z 1/H(z'){\rm d}z'$, and $ R(z) $ denotes the time evolution of the burst rate and has the form

The strain of GWs in GW interferometers can be expressed as [64]

where ψ is the polarization angle, ($ \theta,\phi $) are the angles describing the location of the source in the sky, and $ F_+ $ and $ F_{\times} $ represent the detector antenna pattern functions of the ET [62], which are respectively given by

Because the three interferometers are arranged in an equilateral triangle, the other two antenna pattern functions can be expressed as $ F_{+,\times}^{(2)}(\theta,\phi,\psi)=F_{+,\times}^{(1)}(\theta,\phi+2\pi/3,\psi) $ and $ F_{+,\times}^{(3)}(\theta,\phi,\psi)=F_{+,\times}^{(1)}(\theta,\phi+4\pi/3,\psi) $, respectively.

By applying the stationary phase approximation, the Fourier transform $ \mathcal{H}(f) $ of the strain $ h(t) $ is [62, 65]

with $ \Psi(f) $ being a phase, and the Fourier amplitude is given by

where $ \mathcal{M}_c $ is the chirp mass defined by $ \mathcal{M}_c=M\eta^{3/5} $, $ M=m_1+m_2 $ is the total mass of the coalescence of the binary system with component masses $ m_1 $ and $ m_2 $, and $ \eta=m_1 m_2/M^2 $ is the symmetric mass ratio. In the observer frame, the observed chirp mass is related to the physical chirp mass by a redshift factor, i.e., $ \mathcal{M}_{c,\text{obs}} = (1 + z)\mathcal{M}_{c,\text{phys}} $. Moreover, $ \mathcal{M}_c $ in Eq. (34) denotes the observed chrip mass, whereas ι represents the angle of inclination of the binary's orbital angular momentum relative to the line of sight. The definitions for the phase $ \Psi(f) $ are available in Refs. [62, 65−67]. Given the expected strong beaming of short gamma-ray bursts (SGRBs), the coincident observations of SGRBs suggest that the binaries are predominantly face-on oriented, meaning $ \iota \simeq 0 $, with the maximum inclination being approximately $ \iota = 20^{\circ} $. Notably, averaging the Fisher matrix over the inclination ι and the polarization ψ, subject to the constraint $ \iota < 20^{\circ} $, is nearly equivalent to assuming $ \iota=0 $ in the simulation [65]. Consequently, we can simplify the simulation of GW sources by setting $ \iota=0 $.

The performance of a GW detector is characterized by the one-sided noise power spectral density $ S_h(f) $ (PSD). We adopt the noise PSD of the ET from Ref. [62]. The combined signal-to-noise ratio for a network of three independent interferometers is given by

where $ \rho^{(i)}=\sqrt{\langle\mathcal{H}^{(i)},\mathcal{H}^{(i)}\rangle} $. The inner product of two functions $ a(t) $ and $ b(t) $ is defined as

where $ \tilde{a}(f) $ and $ \tilde{b}(f) $ are the Fourier transforms of $ a(t) $ and $ b(t) $, respectively. The lower limit in frequency for ET is $ f_{\text{lower}}=1 $ Hz, and the upper limit is given by $ f_{\text{upper}} = 2/(6^{3/2}2\pi M_{\text{obs}}) $, where $ M_{\text{obs}} = (1 + z)M_{\text{phys}} $ is the observed total mass.

The standard Fisher matrix method is employed to estimate the instrumental error in the luminosity distance, assuming that this parameter is uncorrelated with any other GW parameters [62, 65]. This results in

Because $ \mathcal{H}\propto d_L^{-1} $, we have $ \sigma_{d_L}^{\text{inst}}\simeq d_L/\rho $. To account for the effect of inclination ι, where $ 0^{\circ}<\iota<90^{\circ} $, we add a factor of $ 2 $ to the instrumental error, resulting in

An additional error due to gravitational lensing must be considered. For the ET, this error is $ \sigma_{d_L}^{\text{inst}}=0.05 z d_L $. Thus, the total uncertainty on luminosity distance is given by

With this information, we can derive the key parameters for GW events, including z, $ d_L $, and $ \sigma_{d_L} $. Consequently, we can simulate 1000 GW events expected to be detected by the ET during its 10-year observation period [68−77] (see Refs. [78, 79] for recent progress).

To estimate the observational constraints on the free parameters of the EMSG models, we employ the MCMC method. The likelihood function for the mock dataset of GW standard sirens is formulated as follows:

where the datasets $ d_L^{\text{obs}}(z_i) $ denote the luminosity distance of the $ 1000 $ simulated GW events with their associated uncertainties $ \sigma_{d_{L,i}} $, and $ d_L^{\text{th}}(z_i) $ is the theoretical prediction on each coalescing event. The fiducial values are set to $ \Omega_{m}^{\text{fid}}=0.315 $ and $ H_0^{\text{fid}}=70 $ km s$ ^{-1} $ Mpc$ ^{-1} $.

Galaxy surveys are directly sensitive to the combination $ f\sigma_8(z) $, which is a nearly model-independent estimator of the observed growth history of the universe. Consequently, many surveys, such as, 2MRS, VVDS, SDSS, 6dFGS, BOSS, and Vipers, offer measurements for this specific estimator. In recent papers [55, 56], a compilation of the current measurements for $ f\sigma_8(z) $ has been prdovided, and we utilize these measurements for our analysis.

To test the two power-law models described above with the datasets, we use Bayesian inference. The $ \chi^2 $ function for the datasets is defined as

where $ V^i=f\sigma_{8,i}-f\sigma_8(z_i;\Theta) $ with Θ the parameters introduced in the model ($ \Theta=H_0,\Omega_{m},\eta_t,\sigma_8 $ for $ n=1/2 $ case and $ \Theta=H_0,\Omega_{m},\Omega_Q,\sigma_8 $ for $ n=1 $ case). Here, $ f\sigma_8(z_i;\Theta) $ is the theoretical growth function given by the models with a set of parameters Θ at redshift $ z_i $, and $ f\sigma_{8,i} $ corresponds to each of the observed data points. Following Ref. [56], the covariance matrix is constructed as consisting of two parts, i.e., the diagonal and nondiagonal parts. The diagonal part is defined as usual:

whereas the nondiagonal part assesses the potential impact of correlations among different data points. We consider the introduction of randomly selected nondiagonal elements in the covariance matrix while maintaining its symmetry. To be explicit, positive correlations are introduced in nine randomly selected pairs of data points. The positions of the nondiagonal elements are selected randomly, and the magnitude of the randomly selected covariance matrix element, denoted as $ C_{ij} $, is set according to

where $ \sigma_i $ and $ \sigma_j $ represent the published 1σ errors of the respective data points $ i, j $.

The prior is a uniform distribution, where $ \Omega_{m} \in $ [0, 0.5], $ h \in [0.6, 0.8] $, $ \eta_t \in [-1, 1] $, and $ \sigma_8 \in [0, 1] $ for $ n=1/2 $. For $ n=1 $, the intervals are $ \Omega_{m} \in [0, 0.5] $, $ h \in [0.6, 0.8] $, $ \Omega_Q \in [-1, 1] $, and $ \sigma_8 \in [0, 1] $. Moreover, the uniform priors are assumed to be the same for the standard sirens, RSD, and combined data.

In this subsection, we explore the constraints on the power-law models through a statistical analysis involving three categories of data combinations: GW, RSD, and GW+RSD. We employ the MCMC method for a minimum likelihood $ \chi^2 $ fit, and the fitting results are summarized in Table 1. The marginalized probability distribution for each parameter ($ \Omega_{m}, H_0, \eta_t, \sigma_8 $ for $ n=1/2 $; $ \Omega_{m}, H_0, \Omega_Q, \sigma_8 $ for $ n=1 $) along with the marginalized 2D confidence contours are depicted in Figs. 1(a) and 1(b).

Data	$ \Omega_{m} $	$ H_0 $	$ \eta_t $	$ \sigma_8 $
GW+RSD	$ 0.283^{+0.033}_{-0.038} $	$ 67.80\pm 0.81 $	$ 0.008\pm 0.03 $	$ 0.73\pm 0.022 $
GW	$ 0.23^{+0.044}_{-0.056} $	$ 68.18\pm 0.86 $	$ 0.096^{+0.071}_{-0.08} $	$ - $
RSD	$ 0.181^{+0.08}_{-0.1} $	$ - $	$ 0.106^{+0.064}_{-0.12} $	$ 0.691^{+0.027}_{-0.036} $
Data	$ \Omega_{m} $	$ H_0 $	$ \Omega_Q $	$ \sigma_8 $
GW+RSD	$ 0.271\pm 0.019 $	$ 68.06\pm 0.62 $	$ 0.00072\pm 0.0006 $	$ 0.717^{+0.025}_{-0.022} $
GW	$ 0.263\pm 0.027 $	$ 67.98\pm 0.72 $	$ 0.0019^{+0.0012}_{-0.0014} $	$ - $
RSD	$ 0.253^{+0.044}_{-0.037} $	$ - $	$ 0.0008\pm 0.001 $	$ 0.694\pm 0.041 $

Table 1.  Marginalized $ 1\sigma $ uncertainties of the parameters for the two viable models in the EMSG theory. Note that, here, $ H_0 $ is in units of km s$ ^{-1} $ Mpc$ ^{-1} $.

Figure 1.  (color online) Constraints on the EMSG cosmological model with the RSD and GW data. The upper and lower panel respectively correspond to the posterior distributions of $ \Omega_m $, $ H_0 $, $ \eta_t $, and $ \sigma_8 $ for $ n=1/2 $ and $ \Omega_m $, $ H_0 $, $ \Omega_Q $, and $ \sigma_8 $ for $ n=1 $.

Let us begin with $ n=1/2 $. Figure 1 presents the fitting results obtained from the RSD datasets, showcasing the fit results: $ \Omega_{m}=0.181^{+0.08}_{-0.1} $, $ \eta_t=0.106^{+0.064}_{-0.12} $, $ \sigma_8=0.691^{+0.027}_{-0.036} $, and no constraint on $ H_0 $. To facilitate a comprehensive comparison, Fig. 1 also includes the best-fit values of the parameters for the simulated GW data are $H_0=68.18\pm 0.86$ km s$ ^{-1} $ Mpc$ ^{-1} $, $ \Omega_{m}=0.23^{+0.044}_{-0.056} $, and $ \eta_t=0.096^{+0.071}_{-0.08} $, which are nearly determined by the fiducial values as mentioned by many works in the literature [68−72]. However, the standard siren observation will improve the cosmological parameter estimation because the parameter degeneracy observed in the RSD data differs from that in the standard sirens. This distinction becomes even more apparent when examining the joint GW+RSD analysis, combining standard probe datasets in the GW and RSD observations. The fit values for the parameters in this combined analysis are $H_0=67.80\pm 0.81$ km s$ ^{-1} $ Mpc$ ^{-1} $, $ \Omega_{m}=0.283^{+0.033}_{-0.038} $, $ \eta_t=0.008\pm 0.03 $, and $ \sigma_8=0.73\pm 0.022 $. Consequently, the shift in the best-fitted parameters, accompanied by a significantly reduced allowed region, demonstrates the potential of leveraging future GW observations to enhance the precision of model parameters in cosmology.

Table 1 explicitly presents the best-fit values of the parameters along with their $ 1\sigma $ uncertainties derived from three distinct data combinations. Note that in the realm of modified gravity, the best-fit values for the matter density parameter $ \Omega_{m} $ and Hubble parameter $ H_0 $ are remarkably consistent with the current estimates within the ΛCDM cosmology framework, as provided by the Planck data release [80]. Furthermore, the deviation parameter $ \eta_t $ resulting from the analysis of three data combinations (RSD, GW, and RSD+GW) indicates that the ΛCDM model ($ \eta_t=0 $) still falls within a $ 68.3\% $ confidence level. However, note that the best-fit value of $ \eta_t $ appears to be slightly larger than $ 0 $, suggesting a potential that the ΛCDM model might not be the most favored cosmological model according to the current observations.

We now investigate parameter inference capabilities for the standard siren and RSD for $ n=1 $; the fit values of the parameters are obtained as $ \Omega_{m}=0.253^{+0.044}_{-0.037} $, $ \Omega_Q=0.0008\pm 0.001 $, $ \sigma_8=0.694\pm 0.041 $, and no constraint on $ H_0 $ for the RSD dataset. The influential role of standard sirens in breaking degeneracies between model parameters is evident, as illustrated by the clearly defined marginalized $ 1\sigma $ and $ 2\sigma $ contours for each parameter in Fig. 1(b). Fitting the $ n=1 $ model to the combined GW+RSD observations yields parameter values of $ H_0=68.06\pm 0.62 $ km s$ ^{-1} $ Mpc$ ^{-1} $, $ \Omega_{m}=0.271\pm 0.019 $, $ \Omega_Q=0.00072\pm0.0006 $, and $ \sigma_8= 0.717^{+0.025}_{-0.022} $. Table 1 summarizes the best-fit values for the three combined datasets. Notably, the deviation from ΛCDM is also slightly larger than $ 0 $, although it is still encompassed within a 68.3% confidence level, which also implies that ΛCDM may not be the most favored cosmological model according to current and future high-redshift GW+RSD observations.

In Fig. 2, we depict the evolution of $ f\sigma_8 $ with respect to the redshift. The blue, orange, and green lines correspond to the power-law models with $ n=1/2 $, $ n=1 $, and ΛCDM models, respectively. Clearly, the value of $ f\sigma_8 $ for $ n=1/2 $ and $ n=1 $ is almost the same for lower redshift, both of which are less than that of ΛCDM model. This indicates that the power-law models of the EMSG theory predicts the slower growth rate of matter fields. For larger values of the redshift, the evolution of the $ f\sigma_8 $ for $ n=1/2 $ decreases slower than that for $ n=1 $ and ΛCDM models. Thus, more observations are required to favor one of them.

Figure 2.  (color online) Evolution of the $ f\sigma_8 $ as a function of redshift z. Here, the observational data are also shown.

Finally, we must make some critical clarifications. In this paper, the RSD data we used are real observational data, but the GW data are simulated. The RSD data are entirely incapable of constraining the Hubble constant, resulting in a clear parameter degeneracy within the parameter space. The GW data can measure the Hubble constant very well. Therefore, when the RSD data are combined with the GW data, the original degeneracy is broken, and the improvement in parameter constraints is very significant. Because the GW data are simulated, the central values of the cosmological parameters are fixed in the simulation; thus, in our final joint fit results, the significance of the central values is not particularly great, and the most meaningful result is the range of errors. Our work aims to demonstrate is this point: future GW observations, if combined with RSD observations to constrain EMSG theory, can significantly enhance the precision of cosmological parameter constraints.

In this paper, we have developed the theory of the linear cosmological perturbations for the energy momentum squared tensor gravity theory. We focus on the scalar and tensor modes becaue the vector mode decays rapidly with the cosmological expansion. For the tensor perturbation, we show that the propagation speed of gravitational wave is the same as the light, which agrees with the recent multimessenger measurements. Regarding scalar perturbation, our focus has been exclusively on the sub-horizon limit. This approach aims to elucidate the equations governing the linear formation of cosmological structures, essential for assessing these theories against existing cosmological precision data. Notably, these equations have been absent in the existing literature, and our investigation seeks to fill this gap.

Subsequently, by utilizing the simulated data generated by the third-generation GW detector, i.e., the Einstein Telescope, we explore the constraining capabilities of GW events on the EMSG theory. We focus on two specific pow-law models, serving as significant viable models for cosmology. Our findings demonstrate that the sensitivity achievable by the ET detector or a comparable third-generation interferometer is sufficient to enhance the current parameter estimates within the EMSG theory. Moreover, we also explore the constraining capabilities of the growth of matter density perturbation with the RSD datasets. Owing to the non-conservative nature of the theory, the evolution equation for the matter density perturbation undergoes modification within these models. We find that the RSD datasets alone may not offer the constraints comparable to the GW observations. However, when combined with the standard sirens, the precision of the constraint is largely imporved because the degeneracies among parameters are mitigated. Although the ΛCDM model (corresponding to $ \eta_t=0 $ or $ \Omega_Q=0 $) still falls within a $ 68.3\% $ confidence level, the nonvanishing of these two modified gravitational parameters indicates that the ΛCDM model might not be the most favored cosmological model according to the current observations. The plot of $ f\sigma_8 $ demonstrates that, for lower redshifts where the data availability is higher, the pow-law models of the EMSG theory predicts a slower growth rate of matter fields. For larger redshifts, the evolution of $ f\sigma_8 $ for the pow-law models decreases slower than the prediction of the ΛCDM model. Thus, more observations would be required to favor one of them.