The properties of metric perturbations are determined in the context of an expanding Universe governed by a modified theory of gravity with a nonminimal coupling between curvature and matter. We analyse the dynamics of the 6 components of a general helicity decomposition of the metric and stressenergy perturbations consisting of scalar vector and tensor sectors. The tensor polarisations are shown to still propagate luminally in agreement with recent data from gravitational interferometry experiments while their magnitude decays with an additional factor sourced by the nonminimal coupling. We show that the production of these modes is associated with a modified quadrupole formula at leading order. The vector perturbations still exhibit no radiative behaviour although their temporal evolution is found to be modified with spatial dependence remaining unaffected. We establish that the scalar perturbations can no longer be treated as identical. We investigate the scalar sector by writing the modified model as an equivalent twofield scalartensor theory and find the same scalar degrees of freedom as in previous literature. The different sectors are paired with the corresponding polarisation modes which can be observationally measured by their effects on the relative motion of test particles thus providing the possibility of testing the modified theory and constraining its parameters.
['General Relativity and Quantum Cosmology']
Gravitational wave polarisations in nonminimally coupled gravity
{'Gravitational wave polarisations in nonminimally coupled gravity': "∗ †  \nMiguel Barroso Varela and Orfeu Bertolami Departamento de F'ısica e Astronomia, Faculdade de Ciˆencias,  \nUniversidade do Porto, Rua do Campo Alegre s/n, 4169-007 Porto, Portugal  \n(Dated: December 18, 2024)  \nThe properties of metric perturbations are determined in the context of an expanding Universe governed by a modified theory of gravity with a non-minimal coupling between curvature and matter. We analyse the dynamics of the 6 components of a general helicity decomposition of the metric and stress-energy perturbations, consisting of scalar, vector and tensor sectors. The tensor polarisations are shown to still propagate luminally, in agreement with recent data from gravitational interferometry experiments, while their magnitude decays with an additional factor sourced by the nonminimal coupling. We show that the production of these modes is associated with a modified quadrupole formula at leading order. The vector perturbations still exhibit no radiative behaviour, although their temporal evolution is found to be modified, with spatial dependence remaining unaffected. We establish that the scalar perturbations can no longer be treated as identical. We investigate the scalar sector by writing the modified model as an equivalent two-field scalar-tensor theory and find the same scalar degrees of freedom as in previous literature. The different sectors are paired with the corresponding polarisation modes, which can be observationally measured by their effects on the relative motion of test particles, thus providing the possibility of testing the modified theory and constraining its parameters.", 'I. INTRODUCTION': "Since their prediction in the early 1900s, gravitational waves (GWs) were seen as a topic of theoretical interest with no hope of any observational detection, mostly due to their significantly small effect on most objects. This all changed with the first detection of GWs by the LIGO experiment in 2015 [1]. Since their experimental confirmation, GWs are now taken as one of the leading opportunities for probing distant events with gravitational effects far beyond any of those that currently take place in our Solar System [2, 3]. This has led to a surge in interest in their theoretical properties in General Relativity (GR) and modified theories, as current and future experiments may provide the means to measure these and thus shed light on a putative more encompassing theory of gravity [3].  \nAs established in their original prediction, GWs in GR are expected to be present in the form of two massless tensor polarisations that propagate luminally through the vacuum. However, modified gravity theories often lead to the presence of additional polarisations, as well as different propagation properties for the tensor modes [4]. For instance, higher-order effective field theories of gravity have been shown to predict deviations from luminal propagation of the tensor modes in the vicinity of Schwarzschild [5] and Reissner-Nordstrom [6] black holes, along with additional deviations in the cosmological context [7]. Many of these properties have been constrained observationally in Earth-based experiments, which have provided us with an upper bound on the mass of the  \ngraviton [8], along with lower and upper bounds on the speed of these waves [9]. These measurements can be used to tightly constrain parameters of modified theories or to even rule them out, which is a complex task in other contexts due to the sensible nature of gravitational effects at local scales [10, 11]. Although measuring relatively small modifications to GW properties around exotic backgrounds such as black hole spacetimes seems unlikely in the near future, the cosmological large-scale dynamics may provide promising opportunities to analyse modifications to GR in the proximity of our planet [12].  \nBesides analysing alterations to the already predicted tensor modes, the presence of additional polarisations of GWs in the context of modified gravity has been extensively researched in past literature [4, 10, 11, 13, 14]. This follows from the modification of the dynamics of the scalar and vector modes one considers when decomposing the metric perturbations into their different helicity sectors [15, 16]. The investigation of the properties of these additional modes can be achieved by an explicit perturbative approach or by applying the Newman-Penrose tetrad formalism [17], as seen, for example, in Refs. [18, 19]. In a closer parallel to the work presented here, it has been shown that minimally coupled f ( R ) theories predict the existence of scalar gravitational waves [20-22], which served as motivation for deeper research into this topic.  \nIn this work, we consider the properties of metric perturbations in a cosmological expanding background governed by a modified theory of gravity with non-minimal coupling (NMC) of matter and curvature [23]. This advances the work conducted in Ref. [18], where only constant curvature backgrounds were considered as an initial probe for the properties of these waves in the modified regime. This kind of NMC theory has also been exten-  \nsively researched in the context of mimicking dark matter profiles [24, 25], analysing the modified theory with solar system constraints [26-28], sourcing cosmological inflation in the early Universe [29-31] and the creation of large-scale structure [32]. More recently, it has been shown that we may solve the open problem of the 'Hubble tension' [33] by considering the late-time effects of the NMC on the model-dependent evolution of cosmic microwave background data used for indirect measurements of the Hubble parameter, while simultaneously providing an observationally compatible mechanism to source the accelerated expansion of the Universe [34]. This further motivates the analysis of the propagation of NMC gravitational waves in the context of an expanding Universe, as it allows for testing of the mathematical consistency of the theory at these scales, as well as providing the means for additional observational tests of the modifications to GR [10, 12].  \nFor simplicity of the analysis carried out here, whenever considering an expanding Universe we use comoving time η , which we relate to cosmic time as dt = a ( t ) dη , and write the background FLRW metric in Cartesian coordinates as  \nds 2 = a 2 ( η )( -dη 2 + dx 2 + dy 2 + dz 2 ) (1)  \nor more concisely as ¯ g µν = a 2 ( η ) η µν . We define the Hubble parameter by H = ˙ a/a = a ' /a 2 as usual, with dots representing derivatives with respect to cosmic time t and primes representing derivatives with respect to comoving time η . We shall mostly use the 'comoving Hubble parameter' H = a ' /a = aH as this is useful for our analysis, although one can find instances in the literature where H and ˙ H are used for easier interpretation in the context of the expanding Universe. Additionally, we define the full metric as  \ng µν = ¯ g µν + a 2 ( η ) h µν = a 2 ( η )( η µν + h µν ) , (2)  \nwhich will simplify the form of the equations and the determination of the propagation speed of the perturbations in case they exhibit radiative behaviour. As done in this equation, we will write background quantities with a bar as ¯ ρ whenever the distinction is necessary, while often omitting this when a quantity is already multiplied by another that is explicitly of perturbative first-order. We use the ( -, + , + , +) signature and set c = 8 πG = 1 for simplicity.  \nWe have organised this paper as follows. We present the nonminimally coupled model, the relevant field equations and their linearly perturbed form in Section II. The helicity decomposition of the metric and stress-energy tensors, along with their associated polarisations and observable effects on the relative motion of test particles, are described in Section III. This is followed by the analyses of the dynamics of the perturbations in GR and in the nonminimally coupled theory, which are included in Sections IV and V respectively. We present an alternative method to determine the polarisation spectrum  \nbased on an equivalent scalar-tensor theory formulation of the NMC model in Section VI. We conclude the work in Section VII, where we discuss the obtained results along with possible extensions of our considerations.", 'A. Action and field equations': "The nonminimally coupled f ( R ) model can be written in action form as [23]  \nS = ∫ dx 4 √ -g [ 1 2 f 1 ( R ) + [1 + f 2 ( R )] L m ] , (3)  \nwhere f 1 , 2 ( R ) are functions of the scalar curvature R , g is the metric determinant and L m is the Lagrangian density for matter fields [23]. The dynamics of General Relativity can be recovered by setting f 1 = R and f 2 = 0. The inclusion of a cosmological constant can be considered by choosing f 1 = R -2Λ. By varying the action with respect to the metric g µν we obtain the field equations [23]  \n( F 1 +2 F 2 L m ) G µν =(1 + f 2 ) T µν +∆ µν ( F 1 +2 F 2 L m ) + 1 2 g µν ( f 1 -F 1 R -2 F 2 R L m ) , (4)  \nwhere we have defined ∆ µν ≡ ∇ µ ∇ ν -g µν □ and F i ≡ df i /dR . As can be seen from the presented equations, the choice of the Lagrangian density is non-trivial, as it directly affects the resulting dynamics via the nonminimal coupling, while only appearing in the form of the related stress-energy tensor in General Relativity [35, 36]. With this in mind, throughout this work, we follow the arguments discussed in Refs. [36, 37] and take the Lagrangian density to be L m = -ρ , with ρ denoting the energy density.  \nApplying the Bianchi identities ∇ µ G µν = 0 to Eq. (4) leads to the non-conservation law [23]  \n∇ µ T µν = F 2 1 + f 2 ( g µν L m -T µν ) ∇ µ R, (5)  \n̸  \nwhich reduces to the usual stress-energy tensor conservation equation when we set f 2 = 0. This modification of the conservation equation follows directly from the nonminimal coupling of matter and curvature, as seen by its direct dependence on f 2 and its independence from the minimally coupled part of the theory described by f 1 . In this work, as we focus on the effects of the NMC model independently of the minimally coupled f ( R ) model, we set f 1 = R and single out the remaining effects with f 2 = 0.  \nAlthough the form of f 2 will be kept abstract throughout the discussion presented in this work, a useful choice  \nis to consider a general power-law expansion in R . Indeed, writing this as  \n1 + f 2 ( R ) = ∞ ∑ n = -∞ ( R R n ) n , (6)  \nwhere R n can be thought of as setting the scale for which the effects of each term in the series become considerable, we can capture the behaviour of the NMC theory at all orders of the curvature, with positive exponent terms dominating in the early Universe and negative exponent terms coming into play at late-times [34]. This means that we expect the evolution of the perturbations considered in this work to be highly epoch-dependent and for the modifications from GR to be strongest at significantly early and late times of the Universe's expansion.", 'B. Background cosmology in NMC model': "As the behaviour of perturbations will be highly dependent on the background dynamics, it is important to review how the expanding nature of the Universe is modified in the NMC theory. This analysis was originally conducted in Ref. [38], with additional results presented in Refs. [34, 39]. By considering the background metric in comoving coordinates as presented in Eq. (1), we can use the previously discussed field Eqs. (4) to arrive at the modified Friedmann equation  \nH 2 = 1 6 F [ 2(1 + f 2 ) a 2 ρ -6 H F ' -a 2 f 1 + a 2 FR ] . (7)  \nThe spatial components of the field equations yield the modified Raychaudhuri equation  \nH 2 +2 H ' = -1 2 F [ 2 F '' +3 H F ' + a 2 f 1 -a 2 FR +2(1 + f 2 ) a 2 p ] (8)  \nand the modified conservation equation of the perfect fluid stress-energy tensor T µ ν = diag( -ρ, p, p, p ), with ρ  \nand p representing the energy density and pressure respectively, leads to the usual result  \nρ ' +3 H ( ρ + p ) = 0 , (9)  \nwhere the choice of L m = -ρ causes the modifications in Eq. (5) to vanish. This greatly simplifies the background cosmological evolution, as we can consider the usual dynamics of the energy density of all kinds of matter with respect to the scale factor ρ ∝ a -3(1+ ω ) , which depends only on the equation of state parameter ω = p/ρ ( ω = 0 for non-relativistic matter, ω = 1 / 3 for radiation).  \nAt the background level, we thus expect deviations from the standard ΛCDM model, which is already enough to induce modifications on the propagation of GWs, whose properties are dependent on the dynamics of the expansive nature of the Universe. Based on the research conducted in Refs. [24, 34, 39], we know that observational data from galaxy rotation curves and direct measurements of the Hubble and deceleration parameters can be explained by the NMC theory with the inclusion of negative exponent terms in the power series form of f 2 and without need for the inclusion of a cosmological constant. This means that to a reasonable approximation, one can take the background to evolve according to simpler observationally motivated models such as ΛCDM, especially if the numerical evolution of the perturbations proves too computationally expensive to perform in parallel to the simulation of the FLRW background dynamics in the NMC theory.", 'C. Perturbed NMC equations': "Considering the general perturbations to the FLRW background metric, for which we maintain the convention introduced at the start of this paper, and the matter content, the NMC-modified field Eqs. (4) yield  \n( F 1 ,R δR +2 F 2 ,R L m δR +2 F 2 δ L m ) R µν +( F 1 +2 F 2 L m ) δR µν -1 2 g µν F 1 δR -1 2 a 2 h µν f 1 -[ δ ( ∇ µ ∇ ν ) -a 2 h µν □ -g µν δ ( □ ) ] ( F 1 +2 F 2 L m ) -[ ∇ µ ∇ ν -g µν □ ] ( F 1 ,R δR +2 F 2 ,R L m δR +2 F 2 δ L m ) = (1 + f 2 ) δT µν + F 2 T µν δR, (10)  \nwhere we have written F i,R = dF i /dR , with the notation F ' i = F i,R R ' still referring to differentiation with respect to η . Given our previous choice of f 1 = R , from now on we shall set F 1 = 1. The linear perturbation of the covariant derivatives is due to the perturbation of the  \nChristoffel symbols  \nδ Γ ρ µν = 1 2 ¯ g ργ [ ∇ µ ( a 2 h γν ) + ∇ ν ( a 2 h γµ ) -∇ γ ( a 2 h µν )] (11)  \nδ ( ∇ µ ∇ ν ) φ = δ Γ ρ µν ∂ ρ φ (12)  \nfor any scalar quantity φ . For reasons that will become clear later, it is convenient to write the perturbed field equations in the form  \nFδR µ ν + R ν µ δF -1 2 δ µ ν δR +[ δ µ ν □ -∇ µ ν ] δF +[ δ µ ν δ ( □ ) -δ ( ∇ µ ∇ ν )] F = (1 + f 2 ) δ µ ν + F 2 T µ ν δR, (13)  \nwhere F = F 1 + 2 F 2 L m = 1 + 2 F 2 L m and so δF = 2( L m F 2 ,R δR + F 2 δ L m ). Additionally, this shows how the explicit dependence on the Lagrangian density is particularly relevant in the context of perturbation theory, as the matter content itself is perturbed, thus leading to a more complex interplay between the evolution of the small changes in the metric and the stress-energy tensor.", 'A. Metric perturbations': "We aim here to analyse general perturbations to an FLRW background metric. The spatial homogeneity of the system means that all background quantities will be functions of comoving time η only. This is especially important for R = R ( η ), as the background modified theory functions f i and their derivatives will be evaluated at this curvature. Most of the assumptions considered in Ref. [18] are no longer valid, such as ∂ µ F i = ∂ µ L m = 0, as these are now evolving quantities in spacetime. This means that the direct analytic conclusions drawn in that work are no longer possible to reach directly, forcing us to explicitly decompose the gravitational waves into 6 general modes and analyse them individually [16]. Nevertheless, when considering the evolution of the perturbations, we may choose to take the sub-horizon (or highfrequency) approximation, in which the background dynamics may be taken to be static in comparison to the rapidly-varying metric perturbations. We write the general helicity decomposition of h µν as [16]  \nh ηη =2 ϕ h ηi = β i + ∂ i γ h ij = 1 a h TT ij + 1 3 Θ δ ij + ∂ i ε j + ∂ j ε i + ( ∂ i ∂ j -1 3 δ ij ∇ 2 ) λ, (14)  \nwhere we assume that h µν → 0 as we move out to spatial infinity, as one would not expect to feel gravitational effects from an infinitely distant event. We see that we have 6 components from the symmetric h TT ij tensor mode, 3 components for each of the vectors ε i and β i , and 4 scalars ϕ, γ, Θ , λ . This gives a total of 16 modes. Due to the invariance of a ( η ) under spatial rotations, the 4 scalars indeed have helicity 0, the vectors have helicity ± 1 and the tensor has helicity ± 2 [15]. Additionally, we impose the constraints  \n∂ i β i = 0 ∂ i ε i = 0 ∂ i h TT ij = 0 δ ij h TT ij = 0 , (15)  \nwhich give 1/1/3/1 individual equations respectively, totalling 16-6=10 independent variables, consistent with a 4 × 4 symmetric tensor. Note also that we have written spatial indices in equations like ∂ i ε i = 0 with no indication of being raised or lowered, as these are interpreted as 'helicity' indices [15]. These defined quantities are gauge-dependent, allowing for further simplification [15].  \nWe thus define the gauge-invariant quantities  \nΦ ≡ -ϕ + 1 a ∂ η [ a ( γ -1 2 ∂ η λ )] Ψ ≡ 1 6 [ -Θ -∇ 2 λ + H ( γ -1 2 ∂ η λ )] Ξ i ≡ β i -∂ η ε i , (16)  \nwhich yields 6 gauge-invariant functions: 1 for Φ, 1 for Ψ, 3 for Ξ i and 6 for h TT ij , minus 3 for ∂ i h TT ij = 0, minus 1 for δ ij h TT ij = 0 and minus 1 for ∂ i Ξ i = 0. We thus have 6 physical and 4 gauge degrees of freedom [16]. The physical quantities are associated with different modes in the relative motion between two test particles [13] and so introduce detectable modifications to the theory. The remaining gauge freedom can be used to choose the conformal Newtonian gauge [15], for which  \nλ = γ = β i = 0 (17)  \nand so the gauge-invariant quantities can be written as  \nΦ ≡ -ϕ Ψ ≡ -1 6 Θ Ξ i ≡ -∂ η ε i , (18)  \nwhich allows the line element to be written as  \nds 2 = a 2 ( η ) [ -(1 + 2Φ) dη 2 + { (1 -2Ψ) δ ij + ∂ i ε j + ∂ j ε i + 1 a h TT ij } dx i dx j ] . (19)  \nThe scalar sector resembles the Newtonian limit of General Relativity, aiding in the interpretation of the effects of the scalar perturbations Ψ and Φ [15].", 'B. Stress-energy tensor perturbations': "Apart from decomposing the metric perturbations, we may also decompose the stress-energy tensor perturbations into their different helicity components [15]. It is convenient to use δT µ ν , whose components can be written as  \nδT η η = -δρ δT i η = S i + ∂ i S δT i j = δpδ i j +Σ i j , (20)  \nwhere we have defined the anisotropic stress tensor Σ i j . This can in turn be decomposed into scalar, vector and tensor parts as  \nΣ ij = ( ∂ i ∂ j -1 3 δ ij ∇ 2 ) σ +( ∂ i σ j + ∂ j σ i ) + σ TT ij , (21)  \nwith σ TT ij being traceless ∂ i σ TT ij = 0 and transverse in the sense that δ ij σ TT ij = 0. Specifically, in the case of a perfect fluid, such as the one considered in the context of the FLRW metric, the perturbations to the stress-energy tensor can be written as [15]  \nδT η η = -δρ δT i η = -( ρ + p ) v i δT i j = δpδ i j = c 2 s δρδ i j , (22)  \nwhere we have used the definition of the speed of sound in a perfect fluid δp = c 2 s δρ ( c 2 s = 0 for non-relativistic matter and c 2 s = 1 / 3 for radiation). Additionally, we have defined the peculiar velocity of particles due to the perturbations, which we can further decompose into its transverse and longitudinal parts as  \nv i = V i + ∂ i v, (23)  \nwhere V is transverse ( ∂ i V i = 0). Comparing with Eq. (20), we identify S = -( ρ + p ) v , S i = -( ρ + p ) V i and  \nΣ i j = 0. This means that the tensor part of the stressenergy perturbations vanishes for a perfect fluid.  \nDue to T µ ν also being a tensor, we expect it to transform under the same gauge transformations as discussed in the previous section. This motivates us to define gauge-invariant stress-energy perturbations in terms of the already defined quantities. These are [15]  \nδρ ∗ = δρ -3 H ( ρ + p )( v + γ ) v ∗ = v + 1 2 λ ' V i ∗ = V i + β i , (24)  \nwhich in our gauge ( λ = γ = β i = 0) simply give  \nδρ ∗ = δρ -3 H ( ρ + p ) v (25)  \nwith v and V i being gauge-invariant.", 'C. Effects of polarisations on test particles': "The presence of GWs can be inferred from the relative 'stretching/squeezing' of intervals in spacetime. Specifically, this has been used in the LIGO experiment to obtain the first concrete detection of gravitational waves passing through our planet [1]. This observation focused on the more relevant tensorial modes predicted by perturbations to vacuum backgrounds in GR. However, all polarisations discussed here present their respective effects on spacetime intervals and can be associated with different types of effects. To see this, we can analyse perturbations to the geodesic deviation equation due to the metric fluctuations [40]  \nD 2 X µ dτ 2 ≡ u α ∇ α ( u β ∇ β X µ ) = R µ αβσ (1) u α u β X σ , (26)  \nwhere u α is the four-velocity of the observer, X µ is the displacement of two objects travelling along infinitesimally separated geodesics and the superscript '(1)' refers to taking only first-order contributions to the Riemann tensor. By assuming a non-relativistic observer, we can set u α ≈ ( a -1 , 0 , 0 , 0) and thus obtain  \n∂ 2 X i ∂η 2 = -a -2 R i ηjη (1) X j , (27)  \nwhere we have used comoving time η to obtain a simpler relation with our previous calculations. Here we have also taken measurements to be made in the local inertial frame (LIF) of the observer, thus simplifying the lefthand side of the equation to a standard partial derivative [40]. Comoving time is particularly useful for this, as our definition of the full metric in Eq. (2) allowed for perturbations to be directly applied to a 'Minkowski' background, as one would observe in its LIF. This allows us to determine the effects of metric perturbations on slowly moving test particles with the use of the perturbed components of R i ηjη . Aligning the propagation of  \nany waves with the z-axis, we can separate the Riemann tensor into different modes as [13, 41]  \nR i ηjη (1) =   P 4 + P 6 P 5 P 2 P 5 -P 4 + P 6 P 3 P 2 P 3 P 1   ij , (28)  \nwhere we have labelled the 6 independent polarisation modes as P i . In order of increasing i , these describe the longitudinal, vector-x, vector-y, +, × and breathing modes [13], which are shown in Figure 1. By calculating the respective linearised Riemann tensor components for the chosen gauge-invariant metric perturbations, we obtain  \nP 1 = ∂ 2 z Φ+Ψ '' + H (Ψ ' +Φ ' ) P 2 = 1 2 ( H ∂ z Ξ 1 + ∂ z Ξ ' 1 ) P 3 = 1 2 ( H ∂ z Ξ 2 + ∂ z Ξ ' 2 ) P 4 = 1 2 a ( H ' h + + H h ' + -h '' + ) P 5 = 1 2 a ( H ' h × + H h ' × -h '' × ) P 6 = H (Ψ ' +Φ ' ) + Ψ '' , (29)  \nwhere the expanding nature of the Universe is present through the comoving Hubble parameter H and its derivative. The association of the different helicity perturbations with the different polarisation modes serves as yet another reason why the chosen decomposition is useful for observational testing. By considering waves varying much faster than the expansion timescale, we obtain the simpler form  \nP 1 = ∂ 2 z Φ+Ψ '' P 2 = 1 2 ∂ z Ξ ' 1 P 3 = 1 2 ∂ z Ξ ' 2 P 4 = -1 2 a h '' + P 5 = -1 2 a h '' × P 6 = Ψ '' , (30)  \nwhich provides a clear connection between the investigated metric perturbations and the relative motion of test particles. Any detection (or non-detection) of modified behaviour following from the analysis of the NMC theory could provide means to further restrict the theory's parameters or even test for its presence [10, 11, 42, 43].", 'A. Stress-energy tensor perturbations': "As discussed in Section II, the conservation equation in the NMC theory is modified to give  \n∇ µ T µ ν = F 2 1 + f 2 ( δ µ ν L m -T µ ν ) ∂ µ R, (62)  \nwhere we take L m = -ρ . Perturbing this equation to linear order gives  \nδ ( ∇ µ T µ ν ) = F 2 1 + f 2 [( δ µ ν δ L m -δT µ ν ) ∂ µ R + ( δ µ ν L m -T µ ν ) ∂ µ ( δR )] + δ ( F 2 1 + f 2 ) ( δ µ ν L m -T µ ν ) ∂ µ R, (63)  \nwhere the left-hand side is unaltered from GR. Noting that in FLRW all background quantities depend only on η , we see that to zeroth order  \n( δ µ ν L m -T µ ν ) ∂ µ R = ( δ η ν L m -T η ν ) ∂ η R = δ η ν ( -ρ + ρ ) ∂ η R = 0 , (64)  \nwhere we have used our chosen form of L m = -ρ and the ηη component of the diagonal background stress-energy tensor T η η = -ρ . We thus see that the final term in the perturbed conservation equation vanishes for all values of the free index ν . Similarly, by analysing the ν = η  \ncomponent we see that  \nδ ( ∇ µ T µ η ) = F 2 1 + f 2 [ ( δ µ η δ L m -δT µ η ) ∂ µ R +( δ µ η L m -T µ η ) ∂ µ ( δR ) ] = F 2 1 + f 2 [ ( δ L m -δT η η ) ∂ η R +( L m -T η η ) ∂ η ( δR ) ] = 0 , (65)  \nas δT η η = δ L m = -δρ . This shows that the ν = η component of the conservation equation is unaltered by the NMC theory and is thus identical to Eq. (32).  \nFor ν = i we have  \nδ ( ∇ µ T µ i ) = F 2 1 + f 2 [ -δT η i ∂ η R +( δ µ i L m -T µ i ) ∂ µ ( δR )] = F 2 1 + f 2 [ -δT η i ∂ η R +( L m -p ) ∂ i ( δR )] = -F 2 1 + f 2 [ v i ∂ η R + ∂ i ( δR )] ( ρ + p ) ,  \n(66)  \nwhere we again use v i = V i + ∂ i v and so  \nδ ( ∇ µ T µ i ) = -F 2 1 + f 2 ( ρ + p ) V i ∂ η R -∂ i [ F 2 1 + f 2 ( ρ + p )( v∂ η R + δR ) ] . (67)  \nHere we have used the fact that f 2 , ρ and p only depend on η at background level. Similarly to the left-hand side of this equation, we have 2 different terms with different helicities. This means we again have 2 separate equations for each spatial component of the conservation equation. The first of these, associated with the vector sector, is  \n4 H ( ρ + p ) V i +[( ρ + p ) V i ] ' = -F 2 1 + f 2 R ' ( ρ + p ) V i = -ln (1 + f 2 ) ' ( ρ + p ) V i ⇒ V i = V ( s ) i ( z ) (1 + f 2 )( ρ + p ) a 4 , (68)  \nwhere we see that as in GR we can separate the spatial dependence of the vector matter perturbations (here again written with a superscript ( s )) from their temporal evolution, with the latter decaying with an additional factor of (1 + f 2 ) due to the nonminimal coupling.  \nFor the scalar sector we have  \nc 2 s δρ +( ρ + p )Φ + [( ρ + p ) v ] ' = -[ 4 H +ln(1 + f 2 ) ' ] ( ρ + p ) v -F 2 1 + f 2 ( ρ + p ) δR, (69)  \nwhere again we find a contribution of a (1 + f 2 ) factor from the NMC on the evolution of v , along with additional coupling to the scalar metric perturbations via the purely scalar-dependent δR . This is explicitly given by  \nδR = 2 a 2 [ ∇ 2 (2Ψ -Φ) -3 H (3Ψ + Φ) ' -6( H 2 + H ' )Φ -3Ψ '' ] , (70)  \nshowing the added complexity introduced by the NMC theory to the coupling between the scalar quantities δρ , v , Φ and Ψ, which is not present in minimally coupled f ( R ) theories, as these present no alterations to the stressenergy conservation equation [20, 23].", 'B. Scalar perturbations': "The scalar sector describes the effects of the Ψ and Φ perturbations. The non-vanishing components of the perturbed Einstein tensor are  \nδG η η = 6 a 2 H 2 Φ+ 2 a 2 (3 H Ψ ' -∇ 2 Ψ) (40)  \nδG i η = 2 a 2 ( H ∂ i Φ+ ∂ i Ψ ' ) (41)  \nδG i j = 1 a 2 [ δ i j ( 2( H 2 +2 H ' )Φ + 2 H (Φ ' +2Ψ ' ) + 2Ψ '' ) + ( δ i j ∇ 2 -∂ i ∂ j ) (Φ -Ψ) ] (42)  \nwith the respective scalar stress-energy perturbations  \nδT η η = -δρ δT i η = -( ρ + p ) ∂ i v δT i j = δpδ i j = c 2 s δρδ i j . (43)  \nThe ( i, j ) equations can be split into a δ i j term and a ∂ i ∂ j term. These two kinds of terms behave as different helicity components and can thus be separated into independent equations [15]:  \n2( H 2 +2 H ' )Φ + 2 H (Φ ' +2Ψ ' ) + 2Ψ '' = a 2 c 2 s δρ (44)  \nand  \nΦ -Ψ = 0 ⇒ Φ = Ψ , (45)  \nrespectively, where we see that the two scalar perturbations are exactly equal due to the chosen convention in their initial definition in terms of the original perturbation variables. This equality is a direct consequence of the vanishing of σ in the anisotropic stress tensor Σ i j in Eq. (21) when assuming a perfect fluid [15].  \nWe can now analyse the ( i, η ) components of the perturbed field equations with Φ = Ψ. These give  \nH ∂ i Φ+ ∂ i Φ ' = -a 2 2 ( ρ + p ) ∂ i v ⇒H Φ+Φ ' = -a 2 2 ( ρ + p ) v, (46)  \nwhich we can then apply to the ( η, η ) equation  \n2 a 2 ( 3 H 2 Φ+3 H Φ ' -∇ 2 Φ ) = 6 a 2 H ( H Φ+Φ ' ) -2 a 2 ∇ 2 Φ = -3 H ( ρ + p ) v -2 a 2 ∇ 2 Φ = -δρ (47)  \nor in the more illuminating form  \n∇ 2 Φ = a 2 2 ( δρ -3 H ( ρ + p ) v ) = a 2 2 δρ ∗ . (48)  \nThis equation is written fully in terms of gauge-invariant quantities and so is valid in any gauge. It is analogous to the classical Newtonian potential equation ∇ 2 Φ = 4 πGρ [15], here written as ρ/ 2 due to the chosen convention 8 πG = 1 .  \nConsidering the δ i j component of the ( i, j ) equation together with the relation Φ = Ψ, we get  \n2( H 2 +2 H ' )Φ + 6 H Φ ' +2Φ '' = a 2 c 2 s δρ, (49)  \nwhich we combine with the ( i, η ) equation to obtain a gauge-invariant form  \nΦ '' +3(1 + c 2 s ) H Φ ' + ( 2 H ' +(1 + 3 c 2 s ) H 2 ) Φ = a 2 2 c 2 s ( δρ -3 H ( ρ + p ) v ) = a 2 2 c 2 s δρ ∗ . (50)  \nBy observing that the right-hand side of this equation matches that of Eq. (48) with an additional factor of c 2 s , we obtain the full independent 'master' equation for Φ  \nΦ '' -c 2 s ∇ 2 Φ+3(1+ c 2 s ) H Φ ' + ( 2 H ' +(1 + 3 c 2 s ) H 2 ) Φ = 0 , (51)  \n̸  \nwhere we see that Φ obeys a wave-like equation with propagation speed c 2 s . The Φ ' term is a friction term that indicates that the wave has the form (Φ a 3(1+ c 2 s ) / 2 ) instead of simply Φ. The physical perturbation Φ thus propagates as a wave for c 2 s = 0, but also decays with the expansion of the Universe as  \n̸  \nΦ = Ψ ∝ a -3 2 (1+ c 2 s ) ∼ √ ¯ ρ ( c 2 s = 0) . (52)  \n̸  \nDuring the radiation-dominated epoch, we see that the perturbations oscillate while decaying at the same rate as the square root of the background matter density, to which they are inherently connected. After solving for the evolution of Φ we may then determine the behaviour of δρ and v from the previous field and conservation equations [15]. Even though we have found a wave-like equation for the scalar perturbations, this radiative behaviour is only present for c 2 s = 0, meaning that in a matter-dominated epoch, such as the recent past of our Universe, this wave-like evolution of scalar perturbations would not be present, with these instead having a time evolution given by the decoupled Eq. (49) with c 2 s δρ = 0. This equation further simplifies during matter domination, as H 2 + 2 H ' = 0 and c 2 s = 0 during this epoch, leading to the temporal equation for the scalar perturbation Φ '' +3 H Φ ' = 0, which has a solution that is constant in comoving time and one obeying Φ ∝ η -5 ∝ a -5 / 2 [15].  \nAdditionally, by combining the perturbed field and conservation equations for the scalar sector, we may determine the evolution of the matter perturbations δρ in more detail. This is particularly simple in the sub-Hubble limit, with the time scale of perturbations being much  \nsmaller than that of the expanding Universe. In this limit, combining the conservation of the stress-energy tensor with the ( η, η ) component of the field equations gives  \nδρ '' -c 2 s ∇ 2 δρ = 0 , (53)  \nwhich shows how the matter fluctuations also evolve with wave-like behaviour, as one would expect from their intrinsic relation to the scalar metric perturbations. The propagation speed is given by the speed of sound in the corresponding type of perfect fluid, meaning that these matter perturbations propagate at the same speed as the scalar potential Φ (or equivalently Ψ), while also not exhibiting wave-like properties in a matter-dominated background ( c 2 s = 0).  \nAlthough we find the scalar perturbations obey wavelike equations in the presence of matter, the same is no longer true in vacuum ( ρ = δρ = 0). In this case, the metric perturbations have a spatial profile given by the Laplace equation ∇ 2 Φ = 0 and time behaviour given by the unsourced version of Eq. (49). This is expected, as GR only predicts the existence of 2 radiative degrees of freedom [16], which will be discussed in the context of the tensor sector below.", 'C. Vector perturbations': "As we consider Ξ i = Ξ i ( η, z ), the condition ∂ i Ξ i = 0 implies Ξ 3 = 0. As in the GR analysis, the fully spatial components of the perturbed field equations have no vector contribution from the stress-energy terms and thus give us the temporal evolution of the vector sector as  \nF Ξ ' i -2( F 2 ρ ) ' Ξ i +2 H F Ξ i =[ F Ξ i ] ' +2 H F Ξ i = 0 ⇒ Ξ i = R i ( z ) a 2 F , (84)  \nwhich shows that the vector perturbations evolve in comoving time η similarly to how they did in GR, with an additional decay due to the nonminimal coupling term in F .  \nThe ( η, i ) components again determine the spatial evolution of the perturbations  \n1 2 a 2 F ∇ 2 Ξ i = (1 + f 2 ) δT η i + F 2 ¯ T η i δR = (1 + f 2 )( ρ + p ) V i ⇒∇ 2 R i = 2 V ( s ) i ( z ) ,  \n(85)  \nwhere in the final step we have removed the η -dependence from the equations determined above and from the conservation equations. This is precisely the same spatial equation as in GR, leading us to the conclusion that the spatial dependence of the vector perturbations seems to be unaltered by the NMC theory. Additionally, the separation of the η and spatial evolution means that we have no radiative behaviour and thus the NMC theory does not predict wave-like evolution for the vector terms in the metric perturbation. The modification of the temporal evolution of the vector sector could provide an opportunity to observationally test the presence of a nonminimal coupling in the gravitational theory [10, 42]. However, this would involve obtaining observations at a considerable variety of redshifts and throughout a wide variety of systems in order to draw any solid conclusions.", 'D. Tensor perturbations': "The tensor sector of the field equations is made up of the traceless-transverse h TT ij , which we have normalised as h ij = h TT ij /a for simplicity of the final equations. These have the same form as in the standard GR derivation [40], with h TT xx = -h TT yy = h + and h TT yx = h TT xy = h × . The non-zero Einstein tensor perturbations are  \nδG x x = -δG y y = -1 2 a 2 [ ( H 2 + H ' ) h + + □ η h + ] (59)  \nδG x y = δG y x = -1 2 a 2 [ ( H 2 + H ' ) h × + □ η h × ] , (60)  \nwhere we have defined □ η = -∂ 2 η + ∇ 2 . As discussed previously, for a perfect fluid there are no anisotropic perturbations and so the tensor sector has δT µ ν = 0, leading to the final equations  \n□ η h × / + = -( H 2 + H ' ) h × / + = -a 2 R 6 h × / + = -a '' a h × / + , (61)  \nwhich represent a wave-like propagation of the tensor modes with luminal speed c 2 gw = 1. Note that we can think of the right-hand side of the wave equations as an 'effective mass' term with m 2 gw = -( H 2 + H ' ). The sign of this quantity depends on the Ricci scalar R . This curvature is 0 during a radiation-dominated epoch ( a ( η ) ∝ η ) and > 0 in a matter-dominated one ( a ( η ) ∝ η 2 ). Thus it seems that in these stages of the evolution of the Universe we could have a negative graviton mass. However, as discussed in Ref. [7], the presence of such a term merely indicates an effective mass, with a vanishing actual mass of the graviton. This distinction is also relevant when considering the causal structure of these waves, where one can typically take the sub-horizon approximation ∂ 2 z h ∼ k 2 h >> ( a '' /a ) h and retrieve the same luminal speed discussed above [7].", 'V. PERTURBATION DYNAMICS IN NONMINIMALLY COUPLING MODEL': 'For simplicity, when analysing GWs in the modified theory we assume all perturbations to be only functions of η and z , as that should not affect any possible radiative behaviour, which can always be aligned with the  \nz -axis without loss of generality. With this in mind, any second-order spatial derivative terms in scalar quantities like δR choose no preferred direction and thus one can replace ∂ 2 z for ∇ 2 when extending to more general coordinate dependence. For each sector, we find that the corresponding equations are analogous to their GR counterparts, allowing for a simple reconciliation with general dependence on all spatial coordinates.  \nSome analysis of the scalar sector perturbations has been conducted in Ref. [32] in the context of cosmological perturbation theory and the formation of the large-scale structure of the Universe. However, in that work, only approximate behaviour was obtained due to neglecting any time derivatives in comparison with spatial derivatives. While we will sometimes invoke the sub-Hubble regime, we do not make this quasi-static approximation, as we aim to investigate the presence of dynamical gravitational wave-like behaviour in the NMC theory. Another notable difference is the inclusion of vector and tensor sector perturbations, leading to a more complete picture of the evolution of metric fluctuations in this modified theory.', 'B. Tensor perturbations': "The tensor sector of the field equations is relatively simple due to both the zeroth and first-order stressenergy tensor having no tensorial helicity terms. As with the perturbations in GR, the equation for the × term comes from the ( x, y ) component of the field equations, while for the + term an identical equation follows from the ( x, x ) = -( y, y ) components. The evolution of the tensor perturbations is described by  \n(1 -2 F 2 ρ ) □ η h + / × +2( F 2 ρ ) ' h ' + / × = [ 2 H ( F 2 ρ ) ' -( H 2 + H ' )(1 -2 F 2 ρ ) ] h + / × , (71)  \nwhere we note that the □ η operator indicates a luminal propagation speed c 2 tensor = 1. This equation clearly reduces to the GR equivalent when we set f 2 = 0. However, we now also have a friction term due to the nonminimal coupling term ( F 2 ρ ) ' . To simplify this, we may use the fact that  \n(1 -2 F 2 ρ ) h '' -2( F 2 ρ ) ' h ' = Fh '' + F ' h ' = √ F ( √ Fh '' +2 F ' 2 √ F h ' ) = √ F ( √ Fh '' +2 F ' 2 √ F h ' +( √ F ) '' h -( √ F ) '' h ) = √ F ( √ Fh ) '' -√ F ( √ F ) '' h (72)  \nto write the equation as  \n□ η ( ˜ h + / × ) = -( H F ' F + F '' 2 F -F ' 2 4 F 2 +( H 2 + H ' ) ) ˜ h + / × , (73)  \nwhere we have used  \n( √ F ) '' √ F = F '' 2 F -F ' 2 4 F 2 (74)  \nand defined ˜ h + / × = √ Fh + / × . We thus see that the gravitational wave amplitude will be scaled by a η -dependent factor of F -1 / 2 . For late-time modifications with f 2 = ( R/R n ) n ( n < 0), as used in the context of removing the Hubble tension [34] or mimicking dark matter profiles [24], we see that F 2 ρ < 0 and so F = 1 -2 F 2 ρ > 1, meaning that gravitational waves decay faster as F increases into the future and f 2 dominates. If we take the background to be slowly evolving in comparison to the waves and introduce a stress-energy source, we obtain the modified gravitational wave equation  \n□ h TT ij = -16 πG 1 + f 2 F δT ij , (75)  \nwhere we have temporarily restored the gravitational constant. This differs from the GR result by a factor of 1+ f 2 F stemming from the non-minimal coupling. We will later identify this as a parameter (Σ = 1+ f 2 F ) that will also modify the weak lensing predictions from the scalar sector perturbations, which may also be thought of as a rescaling of the gravitational constant ˜ G = G Σ [32].", '1. Generation of gravitational waves in NMC theory': "The production of GWs in GR can be analysed by solving the unmodified version of wave Eq. (75) using the Green's function of the □ operator. For the remainder of this section, we shall simplify the equations by writing δT ij as T ij and consider perturbations varying faster than the background, allowing us to approximate t ≈ η and a ( t ) ≈ 1. By applying the GR conservation equation ∂ µ T µν = 0, one can relate the time derivative of the temporal components of the stress-energy tensor to the spatial derivatives of its spatial components. This leads to the well known quadrupole formula  \nh ij ( t, ⃗x ) = 2 G r I ij ( t r , r ) = 2 G r d 2 dt 2 ∫ d 3 x ' x ' i x ' j T 00 , (76)  \nwhich describes the waves generated by a source at a distance r away and at retarded time t r = t -r . This formula indicates that in GR we do not expect monopoles and dipoles to produce gravitational radiation [40].  \nHowever, in the NMC theory, there are two significant changes. These are the modified dynamical value of the gravitational constant ˜ G and the non-conservation of the stress-energy tensor given in Eq. (5). This means that there are additional steps we must take to derive the form of the waves h ij from the matter content that originates them. We start by applying the Green's function [40]  \nh ij ( t, ⃗x ) = 4 G ∫ dt ' ∫ d 3 x ' δ ( t ' -t r ) | ⃗x -⃗x ' | ( 1 + f 2 F ) T ij ≈ 4 G r ∫ d 3 x ' [( 1 + f 2 F ) T ij ] ( t r , ⃗x ' ) , (77)  \nwhere we have assumed a source at a distance | ⃗x | = r away with ⃗x ≫ ⃗x ' . To simplify this equation, we assume that the background R ( t, ⃗x ) ≈ R ( t ) such that we may move the scaling of the modified constant outside the spatial integral. This means we now have  \nh ij ( t, ⃗x ) = 4 G r ( 1 + f 2 F )∣ ∣ ∣ ∣ t = t r ∫ d 3 x ' T ij ( t r , ⃗x ' ) = 4 ˜ G ( t r ) r ∫ d 3 x ' T ij ( t r , ⃗x ' ) , (78)  \nwhich is what one would expect by considering the rescaled gravitational constant as a 'true' constant. We then integrate by parts and neglect boundary terms as we  \nexpect no contributions at infinity, while using the conservation equation to substitute T ij terms by T 00 terms. If the conservation equation was unaltered with respect to GR, as is the case in minimally coupled f ( R ) theories, then the quadrupole formula for the spin-2 sector would be recovered as  \nh ij ( t, ⃗x ) = 2 ˜ G r I ij ( t r , ⃗x ) , (79)  \nwhere in that case we would have ˜ G = G/F 1 . However, the conservation equation in the NMC theory gives  \n∂ k T kj = -∂ t T 0 j + F 2 1 + f 2 ( g µj L m -T µj ) ∂ µ R = -∂ t T 0 j + F 2 1 + f 2 ( g 0 j L m -T 0 j ) ˙ R = -1 1 + f 2 ∂ t [(1 + f 2 ) T 0 j ] , (80)  \nwhere we have again taken R ≈ R ( t ) and taken the background metric to be flat such that g 0 j = 0. We thus have  \nh ij ( t, ⃗x ) = -4 ˜ G r ∫ d 3 x ' x ' i ∂ k T kj = 4 ˜ G r (1 + f 2 ) ∫ d 3 x ' x ' i ∂ t [ (1 + f 2 ) T 0 j ] . (81)  \nFollowing the same steps we develop this expression as  \nh ij ( t, ⃗x ) = -2 ˜ G r (1 + f 2 ) ∫ d 3 x ' x ' i x ' j ∂ t [ (1 + f 2 ) ∂ k T 0 k ] = 2 ˜ G r (1 + f 2 ) ∫ d 3 x ' x ' i x ' j ∂ t [ (1 + f 2 ) ∂ t T 00 ] = 2 ˜ G r [ I ij + F 2 1 + f 2 ˙ R ˙ I ij ] , (82)  \nwhere we have used the conservation equation again to simplify ∂ k T k 0 as  \n∂ k T k 0 = -∂ t T 00 + F 2 1 + f 2 ( g 00 L m -T 00 ) ˙ R = -∂ t T 00 + F 2 1 + f 2 ( ρ -T 00 ) ˙ R = -∂ t T 00 , (83)  \nas we have chosen the matter Lagrangian to be L m = -ρ = -T 00 . We note that the ˙ R term should only be considered when its dynamics are comparable to those of I ij , as it should otherwise be ignored under the subhorizon approximation. Therefore, under the stated assumptions, we find that the production of gravitational waves is still quadrupolar, obeying a similar equation to the one in GR with a rescaled gravitational constant ˜ G and an additional term due to the modified conservation law, which follows from the nonminimal coupling between matter and curvature. As expected, by setting f 1 = R and f 2 = F 2 = 0, we recover the GR quadrupole formula,  \nthus ensuring a smooth reconciliation with GR in the weak modification regime. However, as discussed later in this work, the NMC theory admits scalar gravitational wave degrees of freedom, as previously discussed in Ref. [18]. This additional degree of freedom predicts gravitational radiation from all multipoles, down to monopoles and dipoles, as in most alternative gravitational theories [44]. Nevertheless, showing that the traceless-transverse modes see no effects from monopoles and dipoles at leading order provides clarity on another fundamental property of what we still expect to be the dominant modes of gravitational radiation in the NMC theory.", 'D. Scalar perturbations': "When considering scalar perturbations in GR, we found the equality Ψ = Φ from the purely spatial compo-  \nnts of the field equations, leading to a greatly simplified analysis of the remaining components. Analysing the same purely spatial components again reveals the possibility of separating these into δ i j and ∂ i ∂ j terms, as expected from the helicity decomposition. For simplicity, we show these in the high-frequency limit, which gives  \nδ i j [ 2(1 + 2 F 2 p )Ψ '' +(1 + 2 F 2 p ) ∇ 2 (Φ -Ψ) -(1 + 2 F 2 p -F ) □ η Ψ+ □ η δF -a 2 (1 + f 2 ) δp ] -∂ i ∂ j [ F (Φ -Ψ) + δF ] = 0 , (86)  \nwhere we have not expanded δF into its full dependence on scalar sector quantities for simplicity and again assumed a perfect fluid with p = c 2 s ρ . By setting F = 1, F 2 = 0 and δF = 0, we recover the high-frequency GR result, as expected. There is a strong dependence on the quantity F p ≡ 1 + 2 F 2 p , which interestingly would be precisely the value of F if we had defined the matter Lagrangian in terms of pressure ( L m = p ). We can again set both of these terms equal to zero independently, with the latter giving  \nΨ -Φ = δ ln F = δF F = -2 F 2 δρ -2 ρF 2 ,R δR F , (87)  \nwhich shows that the NMC, present through F , breaks the equality between the different scalar metric perturbations [32, 45]. This significantly complicates the remaining equations, as we are unable to replace all Ψdependence with Φ-dependence like in GR. Additionally, the ∇∇ δF terms in the perturbed field equations lead to first and second-order derivatives of the density fluctuations δρ , present due to the explicit presence of the matter Lagrangian density in F . However, by rewriting the previous relation as δF = (Ψ -Φ) F and inserting this into the perturbed field Eqs. (13), we may remove most of the complexity from the equations, with stress-energy perturbations now only present in δT µ ν , analogously to GR.  \nThe remaining non-trivial equations are obtained from the ( η, η ) component, the ( i, η ) component and the δ i j term in the purely spatial components of the field equations. The choice of only η and z -dependence means that only the i = 3 component of ( i, η ) is non-trivial, with the remaining ( i, η ) components yielding the same equation in the general case and thus not causing any loss of generality. The δ i j term is the same for all i = j components, as expected. This means we are left with 3 non-trivial field equations, along with 2 stress-energy conservation equations for the scalar sector. This would be consistent with the 5 degrees of freedom from Φ , Ψ , v, δρ and δp . However, as we have assumed a perfect fluid with δp = c 2 s δρ , we are only left with 4 independent degrees of freedom and expect one of the equations to be obtained from the others, as is the case with the equations in GR [15].  \nThe ( i, η ) components of the field equations can again be written as ∂ i ( · · · ) = 0, which together with the condition of vanishing perturbations at infinity leads us to  \nset ( · · · ) = 0. This can be written as  \nF [ H (Φ + Ψ) + (Φ + Ψ) ' ] + F ' (2Φ -Ψ) = -a 2 (1 + f 2 )( ρ + p ) v, (88)  \nwhich reduces to the GR result obtained previously when f 2 = 0, as expected. The ( η, η ) component gives  \n-a 2 (1 + f 2 ) δρ = -F ∇ 2 (Φ + Ψ) + 3 F ' Ψ ' +3 H F ' Φ +3 F ( H ' -H 2 )(Ψ -Φ) +3 H ( F [ H (Φ + Ψ) + (Φ + Ψ) ' ] + F ' (2Φ -Ψ)) , (89)  \nwhere the last term can be simplified with the ( i, η ) equation to give  \nF ∇ 2 (Φ + Ψ) -3 F ' (Ψ ' + H Φ) + 3 F ( H 2 -H ' )(Ψ -Φ) = a 2 (1 + f 2 ) δρ ∗ , (90)  \nwith δρ ∗ = δρ -3 H ( ρ + p ) v still gauge-invariant in the NMC theory. This is the modified form of the 'Poissonlike' Eq. (48), with additional metric terms proportional to F ' and an added factor of (1+ f 2 ) multiplying the matter content of the equation due to the nonminimal coupling. The final term on the left-hand side is proportional to F , which does not vanish in GR ( F = 1). However, in that scenario, we know that Ψ = Φ and the term will still be removed, thus ensuring a reconciliation with GR. In the sub-horizon regime, where we take the perturbations to vary much faster than the rate expansion of the Universe ( ∂ z ∼ k >> H ), this equation simplifies to  \n∇ 2 (Φ WL ) ≡ ∇ 2 (Φ + Ψ) = 1 + f 2 F a 2 δρ ∗ ≡ Σ( a, k ) a 2 δρ ∗ , (91)  \nwhich matches the expression found in Ref. [32]. Note that we have taken all background quantities (such as F ) to evolve with the expansion of the Universe as F ' ∼ H F and thus neglect them when compared to derivatives of perturbative quantities. Similarly to what was done in the tensor sector, we have again defined Σ in analogy with the formalism proposed in Ref. [46], where the importance of this parameter in the context of weak lensing observations was discussed. This definition comes from the relationship between the so-called 'weak lensing potential' Φ WL and the density perturbations. As previously pointed out, this effect can be encapsulated as a rescaling of the gravitational constant [32].  \nAll scalar equations presented so far have had explicit dependence on the perturbed matter content. To remove this, we may recognise that to first-order the stressenergy tensor for a perfect fluid obeys T z z + c 2 s T η η = 0, meaning that doing the same to the corresponding components of the perturbed field equations would remove all matter perturbation content. As before, we choose  \nthe high-frequency limit for simplicity, which gives  \n0 = [ 3(1 + c 2 s ) -F (2 + 3 c 2 s ) ] Ψ '' + [ F (2 + c 2 s ) -2(1 + c 2 s ) ] ∇ 2 Ψ + F Φ '' + [ 1 + c 2 s -F (1 + 2 c 2 s ) ] ∇ 2 Φ . (92)  \nThis is the NMC-modified version of the scalar 'wave' Eq. (51). The main difference between these is the dependence on both scalars Φ and Ψ, which are no longer equivalent. It looks like the combination of two wavelike equations for these scalars. Unlike in GR, the assumption of a matter-dominated Universe with c 2 s = 0 no longer fully removes the ∂ 2 z part of the equation. This is due to the presence of the NMC terms proportional to F 2 ρ . The complexity of the coupled equations means we cannot find concrete analytic solutions even under the high-frequency approximation. Nevertheless, one could consider the weak modification regime, which allows for an additional linear treatment of the NMC effects. However, this would only hold in specific circumstances and is thus left as a discussion in Appendix A.", '1. Detection of scalar polarisations': "Even without an analytical solution for the separate metric perturbations, the analysis conducted here provides a possible method for distinguishing between the standard theory and the f ( R ) modified theory. This follows from the initial discussion on the scalar sector in the NMC theory, where we determined that the presence of the nonminimal coupling breaks the equality Ψ = Φ with the introduction of a δ ln F term. When introducing the scalar polarisation components in the tidal tensor ( P 1 and P 6 ), we found that in the high-frequency (subhorizon) limit these depend on the metric perturbations as  \nP 1 = ∂ 2 z Φ+Ψ '' ∼ -k 2 Φ -ω 2 Ψ P 6 = Ψ '' ∼ -ω 2 Ψ , (93)  \nwhere we note the mixed functional dependence in P 1 . We have also assumed a monochromatic plane wave decomposition of the scalar functions for simplicity in the final step, although one could of course consider more complex Fourier space forms with a mixed frequency spectrum in general. These terms correspond to longitudinal and breathing modes respectively and lead to distinct behaviour from the tensorial + and × modes [16]. In GR, only the massless spin-2 modes are expected to be present, while some modified theories of gravity predict the existence of the aforementioned scalar modes [10, 42]. Particularly, this was already established in the context of minimally and nonminimally coupled f ( R ) theories [18, 20, 22]. This means that the detection of such polarisations would serve as a direct test of modified gravity, although their detection alone would not necessarily distinguish between different modified theories without further determination of their properties. With the above  \nresults, we can go further in our predictions, as we not only find the presence of breathing and scalar modes, but we also obtain a clear distinction between their properties, as we no longer have equivalent temporal and spatial metric perturbations in the scalar sector. If we detect both scalar and longitudinal gravitational wave polarisations and are able to compare their effects, one could theoretically use these observations to isolate the effects of the Ψ perturbation present in both P 1 and P 6 and consequently determine the Φ contribution to P 1 . If these are found to be equivalent, then the nonminimally coupled theory could be ruled out, while evidence of a 'decoherence' of Φ and Ψ would indicate the possible presence of the NMC. Of course, this would not necessarily provide definitive proof of its validity over other alternative theories, but it would nevertheless serve as a remarkable step toward obtaining a more general gravitational theory.  \nDetecting additional polarisations of GWs would require careful synchronization of multiple interferometers, as we would need data from different orientations to be able to disentangle the 6 possible independent polarisations [47, 48]. As discussed in Refs. [42, 49], the responses of ground-based laser interferometers to breathing and longitudinal modes are completely degenerate, meaning that these can not be distinguished by present experiments. This follows from the equal pattern functions obtained for these polarisations in the context of laser interferometer arrays [48]. The emergence of space-based GW detectors, such as LISA [50], leads to the possibility of probing wider ranges of frequencies with minimal interference from the atmosphere [22]. This means that the distinction and independent analysis of the breathing and longitudinal modes, along with the consequent comparison of the two scalar metric perturbations, could become feasible in the coming years [47, 51-53]. Additionally, there are many recent developments on the detection methods for high-frequency gravitational waves (HFGWs) beyond the range of interferometry experiments like LIGO and LISA. These HFGWs can be generated in the very early Universe (see, for instance, Ref. [54]) or high-energy astrophysical events [12, 55], therefore providing a direct connection to extreme conditions in which modified gravity would be indispensable. They would also be in an optimal frequency range for an electromagnetic response [56-58], whose observation may be used to detect and separate the six polarisations, as discussed in Ref. [59], where it was also found that there is a good complementarity between the proposed electromagnetic response data analysis and existing ground-based telescopes.  \nAnother notable avenue for observational tests of modified gravity at the other end of the GW frequency spectrum is the analysis of pulsar timing arrays (PTAs) [48]. This method utilises the highly regular electromagnetic pulsations of distant rotating neutron stars, known as pulsars, to detect low-frequency GWs from the cosmological stochastic background or isolated strong-field events like those caused by super-massive black holes [60]. As  \nlight is emitted from the pulsars and received on Earth, both its start and end points are met with the effects of passing gravitational waves from different points in the Universe, leading to path differences that may be observed as measurable alterations to the otherwise periodic pulsation [61]. Each Earth-Pulsar system then serves as an effective interferometer arm of cosmological scale, allowing for the detection of GWs with frequencies ( ∼ nHz) well below those detectable by LIGO ( ∼ 100Hz) or LISA ( ∼ mHz) [60]. Of course, the large distances to these objects mean that we expect a combination of different gravitational waves to affect the observed result. This stochastic gravitational wave background (SGWB) is a complex composition of different wave polarisations, especially if one considers a general metric theory with 6 possible wave modes. Thankfully, this has been thoroughly researched in past literature, where it has been found that the correlation between data from arrays of many different pulsars can be used to decompose this background into its respective polarisations [61-64]. With the constantly increasing number of pulsars under permanent observation and improvements in the respective data analysis, this method provides yet another chance to individually detect the longitudinal and breathing modes, thus allowing for the separation of the effects from Ψ and Φ [65]. This is especially relevant due to the possible early origin of parts of the SGWB, which could provide a window to epochs with some of the more intense NMC presence. Additionally, the longitudinal nature of the P 1 polarisation leads to an increased effect on the irregularity of the pulsar radiation, with sensitivities up to two orders of magnitude greater than those of tensor modes [61].", 'E. Perturbation dynamics in models with f 2 ( R ) = ( R n /R ) n': "The discussion so far has focused on an agnostic form of f 2 ( R ), chosen to keep all arguments as general as possible. However, it is useful to analyse the dynamics resulting from the modified perturbation equations. For this, we choose an inverse power law form given by f 2 = ( R n /R ) n , where n is a positive integer. This choice ensures the decoupling of curvature and matter for large values of R , which from a cosmological point of view corresponds to the early-time Universe, while allowing for the emergence of modified behaviour near the present. Notably, this kind of model has been shown to mimic the effects of dark matter in galaxy rotation curves [24], yield the observed accelerated expansion of the Universe [39] and resolve the Hubble tension [34]. In the latter case, it was found that the NMC model can approximately recreate the late-time behaviour expansion of the ΛCDM model with H 0 ∼ 73 . 2 km/s/Mpc. This means that for small redshifts we can take the evolution of the background quantities from the results of that work, which greatly simplifies the generation of visualisations for the  \nworked example shown here.  \nIn the remainder of this section, we will focus on models with n = 4 and n = 10, as these were shown to capture different kinds of characteristics of the NMC model in Ref. [34]. The most important function to calculate is F = 1 + 2 F 2 L m , which for these models is given by  \nF = 1 + 2 n ( R n R ) n ρ R ≥ 1 , (94)  \nwhere we can consider suitable values for the NMC 'characteristic scale' constants R n from the results on the Hubble tension presented in Ref. [34], namely R 4 = 4 . 1 × 10 4 and R 10 = 4 . 4 × 10 4 in the chosen unit convention ( c = 1).", '1. Vector perturbations': 'As discussed earlier, the vector perturbations have a temporal evolution given by  \nΞ i ∝ 1 a 2 F = (1 + z r ) 2 F ( z r ) , (95)  \nwhere we have defined the redshift as a ∼ (1 + z r ) -1 in order to avoid confusion with the spatial coordinate z used before. Due to the inverse nature of f 2 , we have F 2 < 0, such that F > 1, thus leading to an overall faster temporal decrease of the perturbations due to the expansion of the Universe, as shown in the left panel in Figure 2. Note that the n = 10 model behaves approximately like GR up until z r ≈ 1 . 75, and presents larger perturbation values than n = 4 for all redshifts, with both models predicting a similar ∼ 20% decrease in the present vector perturbation magnitude when compared to GR.', '2. Tensor perturbations': 'The tensor perturbations not only decay with the expansion of the Universe, but also oscillate with their own characteristic frequency, which is redshifted similarly to what is observed for frequencies of light originating in the far past and detected on Earth ( ω obs = ω i a ( t i ) a ( t obs ) ). However, this decreasing frequency is not affected by the NMC model, instead following from the general expanding features of spacetime in the chosen coordinate system. Similarly, we will factor out the decay over distance ( ∼ 1 /r ), an effect present in both GR and in the NMC model, and regard this as due to the dissipation of the waves as they spread over larger distances. With this in mind, we neglect both of these behaviours and focus on the remaining characteristics of the evolution of the amplitude of these waves. From the previously derived equations, we see that even in GR the amplitude of the tensor perturbations will decay as a -1 ∼ (1 + z r ), which has been factored out from the start in the definition of  \nFIG. 2. Time evolution of the amplitude of vector (top) and tensor (bottom) metric perturbations in GR and in the NMC modified theory with f 2 = ( R n /R ) n . Here the time dependence is shown in terms of the corresponding redshift z r , with any temporal dependence present in both GR and the modified theories being factored out for easier comparison, as pointed out in the y-axis label.  \n<!-- image -->  \nthe perturbed metric. However, in the modified model we obtain oscillatory behaviour for the combination √ Fh , which means that the perturbation h will have an additional temporal decrease given by F -1 / 2 . Putting all of this together, we see that  \nh t ∝ 1 a √ F = 1 + z r √ F ( z r ) , (96)  \nwhere h t denotes the amplitude of the observed tensorial perturbations. The results for the chosen example models are shown in the right panel of Figure 2. Again both models are in agreement with GR for larger redshifts ( z r ⪆ 3), with the n = 10 model maintaining an approximate GR-like behaviour until z r ≈ 1 . 75, after which we see that both models predict a ∼ 10% decrease in the presently observed gravitational wave amplitudes due to the expansion of the Universe in the NMC model.', '3. Scalar perturbations': "The scalar perturbations present more complex equations than the remaining sectors, as seen by the lack of decoupling between the equations for the scalar metric fluctuations Φ and Ψ discussed in the sections above. In fact, even when considering a specific example such as the previously chosen model, we are still unable to obtain analytical solutions for the evolution of these quantities. However, one particular parameter can be predicted in the sub-Hubble limit of Eq. (91), where we defined the 'weak lensing parameter' (Σ) in terms of the relation between the weak lensing potential Φ WL ≡ Φ + Ψ and the density perturbations, similarly to what has been done in past works on cosmological parameters (see Ref. [46] for a review). Distinctly from the vector and tensor sectors, the evolution of this quantity depends not only on F and therefore F 2 ρ , but also on the original NMC function f 2 ( R ). Again applying the generated simulation data from previous work on the Hubble tension [34], we can make predictions on the late-time behaviour of this parameter for different NMC models, as shown in Figure 3.  \nFIG. 3. The evolution of the 'weak lensing parameter' Σ in terms of redshift z r in the NMC modified theory with n = 4 and n = 10. The GR prediction is shown for comparison.  \n<!-- image -->  \nWe find that Σ approaches its GR value (Σ GR = 1) for large redshifts, as expected. In the n = 10 model, this parameter departs from unity later than in the n = 4 model, due to the higher suppression of f 2 with n = 10 for large curvatures. Additionally, Σ consistently decreases almost until the present ( z r = 0), reaching a minimum value of Σ ≈ 0 . 9 before increasing to ∼ 0 . 925, with the turning point occurring around z r = 0 . 2. The n = 4 model is similar in many aspects, such as the initial decrease and consequent turning point. However, there are some notable differences, namely the earlier departure from GR and the earlier turning point around z r = 0 . 7. This latter distinction is of particular importance, as the increase of Σ starts early enough for its value to pass the GR prediction around the same time that it starts increasing  \nin the n = 10 model, reaching an increase of ∼ 10% at present. These are noteworthy results, not only because they provide us with predictions for weak lensing properties of the modified theory which could potentially be observationally tested, but also because they offer a clear distinction between the lower-power n = 4 model and its higher-power n = 10 counterpart, based on their effect on weak lensing at lower redshifts.  \nRecent attempts to constrain the present value of Σ in the context of testing modified gravity theories typically obtain values around 1, with 1 σ regions in the ∼ [0 . 9 , 1 . 2] range at best [66], which places both of the NMC models' predictions within error of observational constraints [66-69]. However, in Ref. [70] it was found from the Dark Energy Survey (DES) year 3 results that the present value of the weak lensing parameter, analysed by its deviation from GR, i.e. Σ 0 ≡ Σ( z r = 0) -1, can be constrained to Σ 0 = 0 . 6 ± 0 . 4 from DES alone and to Σ 0 = 0 . 04 ± 0 . 05 from DES combined with external data, both of which show a greater tension with the prediction from the n = 10 NMC model than with that of the n = 4 model. Similar conclusions were reached in the Planck experiment's 2018 results [71]. This is particularly relevant when considering that these constraints were obtained when taking a parametrisation for Σ of the form  \nΣ par ( z r ) = 1 + Σ 0 Ω Λ ( z r ) Ω Λ ( z r = 0) , (97)  \nwhere Ω Λ represents the dark energy density, which does not need to be constant in general. Although this is not equivalent to the predictions of the NMC theories considered here, taking the deviations of Σ from its GR value to be caused by dark energy is in some ways analogous to what happens in late-time NMC models with dependence on inverse powers of R , which have been shown to yield the same accelerated expansion with no need for the presence of dark energy [34, 39]. Improvements to the accuracy of observational constraints on growth parameters in modified gravity could provide crucial evidence to confirm or rule out the presence of modified theories such as the NMC theory considered in this work.", 'VI. GRAVITATIONAL WAVES IN SCALAR-TENSOR THEORY REPRESENTATION': "An alternative method to determine the spectrum of polarisations expected in an f ( R ) NMC gravity model is to look at an equivalent scalar-tensor formulation of the same theory. In fact, it has been shown in Ref. [72] that a two-field scalar-tensor model provides a suitable choice. This is given in the form of a Jordan-Brans-Dicke theory with a potential  \nS = ∫ d 4 x √ -g [Ω R -V ( γ, Ω) + 2 (1 + f 2 ( γ )) L m ] , (98)  \nwhere γ and Ω are scalar fields and we define  \nV ( γ, Ω) = γ Ω -f 1 ( γ ) = γ (Ω -1) (99)  \nas the potential. Note that we have already assumed the form of f 1 ( x ) = x as our work focuses on the NMC theory. Varying the action with respect to the fields yields γ = R , along with  \nΩ = F 1 ( γ ) + 2 F 2 ( γ ) L m = 1 + 2 F 2 ( γ ) L m , (100)  \nwhile varying with respect to the 'physical' metric gives the field equations  \nΩ ( R µν -1 2 g µν R ) =( ∇ µ ∇ ν -g µν □ )Ω -1 2 g µν V ( γ, Ω) +[1 + f 2 ( γ )] T µν , (101)  \nwhich reduce to the result shown in Section II upon substitution of γ and Ω from the other equations. From now on we will identify γ = R , while keeping Ω as an independent field, as its dependence on the matter Lagrangian L m does not allow us to solve for Ω = Ω( γ ). Taking the trace of the field equations provides a useful relation  \n3 □ Ω+2 V ( R, Ω) -R = 3 □ Ω+ R (Ω -2) = (1+ f 2 ) T, (102)  \nwhich shows that Ω obeys a Klein-Gordon equation. Perturbing the trace equation around some slowly evolving background curvature R 0 ≈ 0 and field Ω 0 , we obtain  \n3 □ δ Ω+(Ω 0 -2) δR = (1 + f 2 ) δT + F 2 TδR, (103)  \nwhile the field equations now read  \nΩ 0 δG µν =( ∇ µ ∇ ν -g µν □ ) δ Ω -1 2 η µν (Ω 0 -1) δR + F 2 T µν δR +(1 + f 2 ) δT µν , (104)  \nwhere we have written the metric perturbations over an approximately flat background as g µν = η µν + h µν . We then choose a gauge in which [18]  \n∂ µ ( h µν -1 2 η µν h -η µν δ Ω Ω 0 ) = 0 , (105)  \nwhich leads to  \nδR µν = ∂ µ ∂ ν δ ˜ Ω -1 2 □ h µν δR = □ δ ˜ Ω -1 2 □ h, (106)  \nwhere we have defined δ ˜ Ω = δ Ω / Ω 0 for simplicity. Note that this is equivalent to the quantity δ (ln F ) = δF/F 0 presented in the previous section, where we found that this is the source of the inequality between the two scalar metric perturbations Φ and Ψ. The field equations are then  \n□ [ Ω 0 h µν +(2Ω 0 F 2 T µν +(Ω 0 -2) η µν ) δ ˜ Ω -( F 2 T µν + 1 2 η µν ) h ] = -2(1 + f 2 ) δT µν , (107)  \nwhile the trace equation gives  \n(Ω 0 F 2 T +2Ω 0 -4) □ δ ˜ Ω+ 1 2 (Ω 0 -2 -F 2 T ) □ h = -(1 + f 2 ) δT, (108)  \nwhere moving background quantities inside of the derivative operators is justified when considering a slowly evolving background, which may be taken to be approximately constant in comparison to the rapidly evolving perturbations. Note that by considering only the tensor sector, i.e. by neglecting the scalar degree of freedom and the trace of the perturbations, we get the same result as in Section V B, with a modified gravitational constant ˜ G = (1 + f 2 ) G/ Ω 0 = (1 + f 2 ) G/F , as expected. Even when assuming the more complex case of the presence of background matter content, it is clear that the NMC theory introduces at least one additional radiative degree of freedom in the form of δ ˜ Ω, which is associated with the non-zero trace of the perturbations as we would expect for a scalar degree of freedom.  \nIndeed, when considering the simpler case where no background matter is present ( ¯ T µν = 0 and Ω 0 = 1), while still allowing for matter perturbations to occur due to the interplay of curvature and matter in the NMC model, the field equations give  \n□ ˜ h µν ≡ □ [ h µν -1 2 η µν h -η µν δ Ω ] = -2(1 + f 2 ) δT µν , (109)  \nwhere the redefined perturbation ˜ h µν is transverse ( ∂ µ ˜ h µν = 0) due to the previously chosen gauge. In the absence of sources, we thus have a wave equation for the quantity ˜ h µν , which is initially sourced by the matter perturbations at the origin of the GW signal, while the trace equation directly shows the relation between δ Ω and h as □ δ Ω = -1 4 □ h . This means that we may solve this equation with the same traceless-transverse solution used in GR. We then invert the definition of ˜ h µν by using h = -˜ h -4 δ Ω to find the physical metric perturbations  \nh µν = ˜ h µν -1 2 η µν ˜ h -η µν δ Ω = ˜ h µν -η µν δ Ω , (110)  \nwhere we have taken ˜ h = 0 due to the traceless property taken for the solution. The physical metric perturbation is then given by  \nh µν =    δ Ω 0 0 0 0 h + -δ Ω h × 0 0 h × -h + -δ Ω 0 0 0 0 -δ Ω    µν , (111)  \nwhich indicates the presence of luminally propagating breathing and longitudinal polarisation modes induced by the scalar perturbation δ Ω along the standard + and × polarisations present in GR [16, 44]. Note that in terms of the scalar perturbations discussed in Section V, this would mean Ψ = -Φ = δ Ω / 2, leading to the exact relation we found from using a gauge-invariant formalism  \nΨ -Φ = δ Ω = δF . We also see that this perturbation obeys h = -4Ω, which agrees with the trace equation in vacuum.  \nHowever, it is important to keep in mind that the NMC contribution to δ Ω follows from two distinct terms, namely F 2 δ L m and F 2 ,R L m δR , meaning that the existence of these waves would depend on the presence of background matter to allow for their propagation. When such waves pass through cosmic voids ( L m = δ L m = 0), as will inevitably be the case in a Universe where dark energy is a purely gravitational phenomenon, both of these contributions are negligible or even zero. The specific case of NMC GW perturbations over dark energy-like fluids has been analysed in Ref. [18], but we do not consider these in this work as the NMC theory has no need for dark energy to generate the accelerated expansion of the Universe [34, 39]. Of course, the prior reasoning only applies if we consider a purely NMC theory with no modification to the minimal sector through f 1 , as minimally coupled f ( R ) theories introduce scalar gravitational radiation polarisations even in vacuum [44], as these do not vanish by taking L m = δ L m = 0. This means that, although the non-detection of scalar polarisations in future experiments could rule out the validity of minimally coupled f ( R ) theories, the same cannot be said about their NMC counterparts, as to leading order the scalar polarisations may have been 'washed out' by the passage of the GWs through cosmic voids.", 'VII. CONCLUSIONS': "In this work, we have analysed the properties of metric and matter perturbations in the context of an expanding Universe in a nonminimally coupled theory of gravity. These same perturbations have been thoroughly studied in the context of GR[15, 40], establishing the evolution of small deviations in homogeneous FLRW spacetimes and matter backgrounds and their role in the formation of large-scale structure in the Universe. The latter studies allow for the prediction of the propagation properties of (tensorial) gravitational waves, which radiate at luminal speeds and decay with the expansion of the Universe, as one would expect [7, 15]. GR also predicts the presence of 2 scalar and 2 vectorial degrees of freedom. While the latter do not exhibit any radiative behaviour, evolving independently with respect to time and spatial coordinates, the former behave as one equivalent term with different epoch-dependent properties [15]. The corresponding scalar sector term obeys a wave-like equation during the radiation-dominated epoch of the Universe, propagating at the speed of sound in that medium ( c 2 s = 1 / 3). However, during matter domination ( c 2 s = 0) the time and spatial dependencies separate, leading to a decaying evolution with no oscillatory behaviour. This means that the only radiative behaviour one would expect to detect in a GR-ruled Universe would come from the usual tensorial degrees of freedom, typically separated into + and  \n× modes [40].  \nWhen considering the corresponding field equations in the NMC theory, several alterations arise [18]. Although the tensorial modes still propagate luminally, their effective mass is altered and the gravitational waves decay with an additional factor of F -1 / 2 while oscillating. These changes are in line with the analogous result presented in Ref. [18], where the additional time dependence was not obtained due to the assumption of simple constant curvature backgrounds. In late-time conditions, we expect negative exponent terms in a full power series expansion of f 2 to dominate [34, 39], thus leading to a faster decrease in the wave magnitude as it propagates. We then analysed the production of these waves, which under our assumptions maintain their quadrupolar nature at leading order, while monopole or dipole sources would only induce the scalar polarisations [44].  \nThe evolution of vector perturbations is mostly unaltered, with spatial and temporal dependencies separating as in GR, except for the latter being modified by a factor of F -1 , which accelerates their attenuation over time. Even though we cannot classify this sector as containing gravitational waves, any effects stemming from these kinds of metric fluctuations would be affected by the nonminimally coupled theory, thus providing means for future testing of the presence of such modifications to GR [10, 42].  \nThe scalar sector has the most non-trivial changes, mostly due to its inherent connection with the scalar density perturbations, which directly enter the field and conservation equations via the nonminimal coupling and due to our choice of Lagrangian density L m = -ρ [23, 35]. This is particularly relevant when considering the relation between the scalar metric perturbations Ψ and Φ, which are no longer simply equivalent as in GR, but rather related by a more complex coupled equation sourced by the variation in F = F 1 + 2 F 2 L m [32]. This quantity depends on the matter density and the Ricci scalar, which are affected by the perturbations to the stress-energy content and metric components respectively. However, the same equation can be used to remove all density perturbation derivative terms from the field equations, thus allowing for a better comparison between the GR and NMC dynamics. We derived the modified 'Poisson-like' equation, where the Laplacian term is now related to the gauge-invariant density perturbations via an 'effective gravitational constant' ˜ G = G (1+ f 2 ) /F which is unique to the NMC theory [32]. We analysed the effect of this modification on the weak lensing parameter (Σ) in examples where f 2 = ( R n /R ) n , with the n = 4 , 10 models chosen from their previously determined behaviour in Ref. [34]. We found our predictions to be in good agreement with recent observational constraints.  \nA particularity of the modified theory considered in this work is the aforementioned inequality of the scalar metric fluctuations Φ and Ψ, as initially introduced in Ref. [32] in the context of large-scale structure formation. Although this adds considerable complexity to the ana-  \ntical determination of properties of the scalar modes, it provides us with a general prediction which could be used in tests of the validity of the NMC theory. We have outlined how this hinges on the detection and distinction of the associated breathing and longitudinal modes, which depend on future space-based gravitational wave interferometry experiments, as the two scalar polarisations produce degenerate responses in ground-based laser interferometer experiments [42, 49]. With the predicted creation of future experiments such as LISA, the status of modified theories like the one addressed in this work could thus be put to the test [42]. Recent interest in the detection of lower and higher frequency signals of the SGWB spectrum could also bring about additional tests of modified gravity theories. Among some of the possibilities discussed in this work are the analysis of PTA data ( ∼ nHz) [60] and the proposed adaptation of microwave cavity experiments ( ∼ MHz) [56]. This is a fundamental step in the analysis of the NMC f ( R ) theory, which has already been shown to fit several observational effects (see, for example, Refs. [24, 34, 39]), while only now providing predictions that can be tested without being led by previous results, as any consistent theory must do.  \nWe then analysed an alternative method to determine the polarisation spectrum in the NMC theory by making use of the representation of this model as a two-field scalar-tensor theory [72]. We find that the scalar field associated with the nonminimal coupling induces an additional degree of freedom that mixes curvature and matter perturbations. The scalar degree of freedom is associated with the trace of the metric perturbations and presents itself in the form of breathing and longitudinal polarisation modes propagating luminally, in agreement with what has been found in f ( R ) theories [44] and in other work on the NMC theory [18]. However, we discussed the particularities of the propagation of these scalar modes through empty regions such as cosmic voids [73], in which the absence of background matter to perturb would remove any NMC effects, in particular the scalar polarisations. This only applies to pure NMC models with no minimal f ( R ) component, meaning that a future lack of detection of scalar GW polarisations from distant sources could provide insight into an experimental distinction between mixed (MC+NMC) and pure nonminimally coupled f ( R ) theories.  \nExtensions of the work described here are varied, with perhaps the clearest direction being performing the analysis of the scalar sector under a general regime. Particularly, research on the evolution of the density perturbations has been conducted in Ref. [32], where the sub-horizon regime was assumed for simplification of the coupled equations. Once a clearer picture of the evolution of scalar quantities is obtained, a logical extension would be to apply those results, together with those presented here, to the prediction of observable primordial imprints of the different perturbation sectors in cosmological data, such as the CMB power spectrum [74]. Additionally, applying the Newman-Penrose tetrad formalism to the  \nperturbed NMC theory in a cosmological context could provide further conclusions on the predicted polarisations and their properties, as done in Refs. [18, 19].", 'ACKNOWLEDGMENTS': 'The work of one of us (O.B.) is partially supported by FCT (Funda¸c˜ao para a Ciˆencia e Tecnologia, Portugal) through the project CERN/FIS-PAR/0027/2021, with DOI identifier 10.54499/CERN/FIS-PAR/0027/2021.', 'Appendix A: Scalar perturbations in the weak NMC regime': "Decoupling the scalar field equations analytically in an abstract NMC theory with an unspecified f 2 ( R ) is extremely convoluted. Employing a numerical method could be necessary to extract more quantitative results in general scenarios [15], which with the results of this work would simply consist of evolving the coupled differential equations with a robust numerical integrator with sensible initial conditions and a set chosen NMC theory parameters. With this in mind, some results on a concrete example were presented in Section V E. Nevertheless, some additional conclusions on the scalar sector can be reached in the weak regime, in which we take all modifications to be treated perturbatively to linear order in a separate expansion to that of the metric perturbations. This assumption is not expected to hold under general conditions, as can be seen by considering the hypothetical complete power expansion form of the nonminimal coupling function f 2 discussed in Section II. The behaviour of negative and positive exponents of the curvature scalar is such that we expect the strength of the NMC to be dominant in spacetime regions with significantly small (late-time Universe) [34, 39] or large (earlytime Universe) [29-31] curvatures respectively. However, intermediary situations, such as those expected around the creation of the cosmic microwave background, could provide regions in which the NMC might be treated in the weak regime, as discussed in Ref. [34]. The following perturbative analysis, although not general, is therefore still a relevant example to consider when describing some fundamental points of the evolution of the Universe.  \nIn order to simplify the distinction between expansion orders, we introduce a new factor χ multiplying f 2 . As discussed in Refs. [5, 6], this method allows us to analyse any modifications to the equations under the light of their  \nGRcounterpart, as they are already first-order terms and thus can be evaluated with the zeroth-order equations in mind. This is particularly useful, as the scalar perturbations Ψ and Φ are equivalent in GR. More specifically, in the perturbative regime we may write  \nΨ -Φ = δF F ≈ δF = -2 χ ( ρF 2 ,R δR + F 2 δρ ) = O ( χ ) , (A1)  \nwhich may be introduced into the field Eqs. (13). To decouple the equations, we need to remove the δρ dependence, which can be done by using the ( η, η ) component of the GR perturbed field equations, as the δρ term above is already O ( χ ). By taking a Fourier wave form for each scalar sector function with ∇ 2 ∼ -k 2 and the explicit expression for δR given previously, we can rearrange the equation above as  \nΨ = Φ + χq (Φ , Φ ' , Φ '' ) (A2)  \nor equivalently  \nΦ = Ψ -χq (Ψ , Ψ ' , Ψ '' ) , (A3)  \nwhere we have introduced the function q (Φ , Φ ' , Φ '' ) for simplicity. This can be applied to the coupled master Eq. (92) by writing Ψ ' = Φ ' + χ ˜ q (Φ , Φ ' , Φ '' , Φ (3) ) and equivalently for Ψ '' , hence removing all Ψ dependence, while now obtaining fourth-order derivative terms in the decoupled master equation.  \nHowever, any third or fourth-order derivatives are necessarily O ( χ ) as the zeroth-order master equation only included derivatives of up to second order. We may thus apply Eq.(51) to lower any higher-order derivatives and write the master equation in the form [7]  \n∂ 2 η Φ+ β 1 ( η ) ∂ η Φ+ β 0 ( η, k 2 )Φ + α 0 ( η )Φ = 0 , (A4)  \nfrom which we can read off the propagation speed c 2 Φ = β 0 /k 2 , with equivalent reasoning leading to analogous equations for Ψ. This is clear when disregarding the friction term, which simply leads to a modified decay of the wave, after which one may take a wave-like time dependence of the form e iωη and solve the equation in the largek limit with ω 2 = β 0 + α 0 . This implies the previously asserted propagation speed c 2 Φ = lim k →∞ ω 2 k 2 = β 0 /k 2 , as the α 0 /k 2 term is negligible in this limit.  \nAs expected, these speeds are of the form c 2 Φ / Ψ = c 2 s + χ ∆ c 2 Φ / Ψ , where these are given by  \n̸  \n- χ\n- c 2 Φ = c 2 s -a 2 [ -2 H 2 ρF 2 ,R +2 ( H 2 + H ' ) ρF 2 ,R -6 H 2 c 2 s ρF 2 ,R -12 ( H 2 + H ' ) c 2 s ρF 2 ,R +36 H 2 c 4 s ρF 2 ,R -18 ( H 2 + H ' ) c 4 s ρF 2 ,R -2 H 2 F 2 +2 ( H 2 + H ' ) F 2 -18 H 2 c 2 s F 2 +6 ( H 2 + H ' ) c 2 s F 2 -a 2 ρF 2 +2 a 2 c 2 s ρF 2 +3 a 2 c 4 s ρF 2 + +2 H F 2 ,R ρ ' -18 H c 4 s F 2 ,R ρ ' +2 H ρF ' 2 ,R -18 H c 4 s ρF ' 2 ,R -4 ρ ' F ' 2 ,R +12 c 2 s ρ ' F ' 2 ,R +2 H F 2 ' +6 H c 2 s F ' 2 -2 F 2 ,R ρ '' +6 c 2 s F 2 ,R ρ '' -2 ρF '' 2 ,R +6 c 2 s ρF '' 2 ,R -2 F '' 2 ] (A5)\n- c 2 Ψ = c 2 s -χ a 2 [ -6 H 2 F 2 +2 ( H 2 + H ' ) F 2 +42 H 2 c 2 s F 2 -18 ( H 2 + H ' ) c 2 s F 2 -a 2 ρF 2\n- +2 a 2 c 2 s ρF 2 +3 a 2 c 4 s ρF 2 -6 H 2 ρF 2 ,R +2 ( H 2 + H ' ) ρF 2 ,R +54 H 2 c 2 s ρF 2 ,R\n- -12 ( H 2 + H ' ) c 2 s ρF 2 ,R -108 H 2 c 4 s ρF 2 ,R +54 ( H 2 + H ' ) c 4 s ρF 2 ,R (A6)\n- -2 H F 2 ,R ρ ' +18 H c 4 s F 2 ,R ρ ' -2 H F ' 2 -6 H c 2 s F ' 2 -2 H ρF ' 2 ,R\n- +18 H c 4 s ρF ' 2 ,R +4 ρ ' F ' 2 ,R -12 c 2 s ρ ' F ' 2 ,R +2 F 2 ,R ρ ''\n- -6 c 2 s F 2 ,R ρ '' +2 F '' 2 +2 ρF 2 ,R '' -6 c 2 s ρF '' 2 ,R ] .  \nNot only do these speeds match at zeroth order, but they still depend on the speed of sound from the perfect fluid considered for the FLRW background. Although in an RD Universe ( c 2 s = 1 / 3) the weak NMC corrections would be relatively small in comparison to the GR speed, this is no longer the case in the more recent MD Universe ( c 2 s = 0), where the radiative behaviour would be fully due to the NMC corrections, even when considering the weak regime. This means that we could in principle use the detection of this behaviour of scalar polarisations of  \nGWsvia their effects on test particles to test the presence of a nonminimal coupling in the gravitational theory. We should note that we have ignored the possibility of a dark energy-dominated (ΛD) Universe, as the source of the currently observed accelerated expansion has been shown to be accounted for in the NMC theory with no need for the presence of a cosmological constant [34, 39]. The difference between the scalar perturbation speeds in the NMC theory ( c 2 Ψ = c 2 Φ ) could provide additional avenues of testing [14].  \n- [1] B. P. Abbott et al. (LIGO Scientific, Virgo), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116 , 061102 (2016), arXiv:1602.03837 [gr-qc].\n- [2] M. Sakellariadou, Gravitational Waves: The Theorist's Swiss Knife, Universe 8 , 132 (2022), arXiv:2202.00735 [astro-ph.CO].\n- [3] L. Xu, Gravitational Waves: A Test for Modified Gravity, Phys. Rev. D 91 , 103520 (2015), arXiv:1410.6977 [astroph.CO].\n- [4] C. Bogdanos, S. Capozziello, M. De Laurentis, and S. Nesseris, Massive, massless and ghost modes of gravitational waves from higher-order gravity, Astropart. Phys. 34 , 236 (2010), arXiv:0911.3094 [gr-qc].\n- [5] C. de Rham, J. Francfort, and J. Zhang, Black Hole Gravitational Waves in the Effective Field Theory of Gravity, Phys. Rev. D 102 , 024079 (2020), arXiv:2005.13923 [hepth].\n- [6] M. Barroso Varela and H. Rauch, Gravitational waves on charged black hole backgrounds in modified gravity, Gen. Rel. Grav. 56 , 16 (2024), arXiv:2311.07376 [gr-qc].\n- [7] C. de Rham and A. J. Tolley, Speed of gravity, Phys. Rev. D 101 , 063518 (2020), arXiv:1909.00881 [hep-th].\n- [8] C. de Rham, J. T. Deskins, A. J. Tolley, and S.-Y. Zhou, Graviton Mass Bounds, Rev. Mod. Phys. 89 , 025004  \n(2017), arXiv:1606.08462 [astro-ph.CO].  \n- [9] N. Cornish, D. Blas, and G. Nardini, Bounding the speed of gravity with gravitational wave observations, Phys. Rev. Lett. 119 , 161102 (2017), arXiv:1707.06101 [gr-qc].\n- [10] B. P. Abbott et al. (LIGO Scientific, Virgo), Search for Tensor, Vector, and Scalar Polarizations in the Stochastic Gravitational-Wave Background, Phys. Rev. Lett. 120 , 201102 (2018), arXiv:1802.10194 [gr-qc].\n- [11] M. Isi, A. J. Weinstein, C. Mead, and M. Pitkin, Detecting Beyond-Einstein Polarizations of Continuous Gravitational Waves, Phys. Rev. D 91 , 082002 (2015), arXiv:1502.00333 [gr-qc].\n- [12] S. Mastrogiovanni, D. A. Steer, and M. Barsuglia, Probing modified gravity theories and cosmology using gravitational-waves and associated electromagnetic counterparts, Phys. Rev. D 102 , 044009 (2020).\n- [13] Y.-Q. Dong, Y.-Q. Liua, and Y.-X. Liu, Polarization modes of gravitational waves in general modified gravity: General metric theory and general scalar-tensor theory, Phys. Rev. D 109 , 044013 (2024), arXiv:2310.11336 [gr-qc].\n- [14] K. Schumacher, N. Yunes, and K. Yagi, Gravitational wave polarizations with different propagation speeds, Phys. Rev. D 108 , 104038 (2023).\n- [15] M. Maggiore, Gravitational Waves. Vol. 2: Astrophysics\n- and Cosmology (Oxford University Press, 2018).\n- [16] E. E. Flanagan and S. A. Hughes, The Basics of gravitational wave theory, New J. Phys. 7 , 204 (2005), arXiv:grqc/0501041.\n- [17] E. Newman and R. Penrose, An Approach to gravitational radiation by a method of spin coefficients, J. Math. Phys. 3 , 566 (1962).\n- [18] O. Bertolami, C. Gomes, and F. S. N. Lobo, Gravitational waves in theories with a non-minimal curvaturematter coupling, Eur. Phys. J. C 78 , 303 (2018), arXiv:1706.06826 [gr-qc].\n- [19] M. E. S. Alves, O. D. Miranda, and J. C. N. de Araujo, Probing the f(R) formalism through gravitational wave polarizations, Phys. Lett. B 679 , 401 (2009), arXiv:0908.0861 [gr-qc].\n- [20] D. Liang, Y. Gong, S. Hou, and Y. Liu, Polarizations of gravitational waves in f ( R ) gravity, Phys. Rev. D 95 , 104034 (2017), arXiv:1701.05998 [gr-qc].\n- [21] H. Rizwana Kausar, L. Philippoz, and P. Jetzer, Gravitational wave polarization modes in f ( R ) theories, in 14th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Astrophysics, and Relativistic Field Theories , Vol. 2 (2017) pp. 1220-1226.\n- [22] T. Katsuragawa, T. Nakamura, T. Ikeda, and S. Capozziello, Gravitational Waves in F ( R ) Gravity: Scalar Waves and the Chameleon Mechanism, Phys. Rev. D 99 , 124050 (2019), arXiv:1902.02494 [gr-qc].\n- [23] O. Bertolami, C. G. Boehmer, T. Harko, and F. S. N. Lobo, Extra force in f(R) modified theories of gravity, Phys. Rev. D 75 , 104016 (2007), arXiv:0704.1733 [gr-qc].\n- [24] O. Bertolami, P. Fraz˜ao, and J. P'aramos, Mimicking dark matter in galaxy clusters through a non-minimal gravitational coupling with matter, Phys. Rev. D 86 , 044034 (2012), arXiv:1111.3167 [gr-qc].\n- [25] O. Bertolami and J. P'aramos, Mimicking dark matter through a non-minimal gravitational coupling with matter, JCAP 03 , 009, arXiv:0906.4757 [astro-ph.GA].\n- [26] R. March, O. Bertolami, J. P'aramos, and S. Dell'Agnello, Nonminimally coupled curvature-matter gravity models and Solar System constraints, in 15th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Astrophysics, and Relativistic Field Theories (2019) arXiv:1903.07059 [gr-qc].\n- [27] R. March, J. P'aramos, O. Bertolami, and S. Dell'Agnello, Nonminimally coupled gravity and planetary motion, Frascati Phys. Ser. 64 , 33 (2017).\n- [28] R. March, O. Bertolami, M. Muccino, and S. Dell'Agnello, Equivalence principle violation in nonminimally coupled gravity and constraints from lunar laser ranging, Phys. Rev. D 109 , 124013 (2024), arXiv:2312.14618 [gr-qc].\n- [29] O. Bertolami, V. Bessa, and J. P'aramos, Inflation with a massive vector field nonminimally coupled to gravity, Phys. Rev. D 93 , 064002 (2016), arXiv:1511.03520 [grqc].\n- [30] O. Bertolami, P. Fraz˜ao, and J. P'aramos, Reheating via a generalized non-minimal coupling of curvature to matter, Phys. Rev. D 83 , 044010 (2011), arXiv:1010.2698 [gr-qc].\n- [31] C. Gomes, J. a. G. Rosa, and O. Bertolami, Inflation in non-minimal matter-curvature coupling theories, JCAP 06 , 021, arXiv:1611.02124 [gr-qc].\n- [32] O. Bertolami, P. Fraz˜ao, and J. P'aramos, Cosmological perturbations in theories with non-minimal cou-  \n- ing between curvature and matter, JCAP 05 , 029, arXiv:1303.3215 [gr-qc].\n- [33] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, In the realm of the Hubble tension-a review of solutions, Class. Quant. Grav. 38 , 153001 (2021), arXiv:2103.01183 [astro-ph.CO].\n- [34] M. Barroso Varela and O. Bertolami, Hubble tension in a nonminimally coupled curvature-matter gravity model, JCAP 06 , 025, arXiv:2403.11683 [gr-qc].\n- [35] O. Bertolami and J. P'aramos, Homogeneous spherically symmetric bodies with a non-minimal coupling between curvature and matter: the choice of the Lagrangian density for matter, Gen. Rel. Grav. 47 , 1835 (2015), arXiv:1306.1177 [gr-qc].\n- [36] O. Bertolami, F. S. N. Lobo, and J. P'aramos, Nonminimal coupling of perfect fluids to curvature, Phys. Rev. D 78 , 064036 (2008).\n- [37] J. D. Brown, Action functionals for relativistic perfect fluids, Classical and Quantum Gravity 10 , 1579 (1993).\n- [38] O. Bertolami and J. P'aramos, Modified Friedmann Equation from Nonminimally Coupled Theories of Gravity, Phys. Rev. D 89 , 044012 (2014), arXiv:1311.5615 [gr-qc].\n- [39] O. Bertolami, P. Fraz˜ao, and J. P'aramos, Accelerated expansion from a non-minimal gravitational coupling to matter, Phys. Rev. D 81 , 104046 (2010), arXiv:1003.0850 [gr-qc].\n- [40] R. M. Wald, General Relativity (Chicago Univ. Pr., Chicago, USA, 1984).\n- [41] D. M. Eardley, D. L. Lee, and A. P. Lightman, Gravitational-wave observations as a tool for testing relativistic gravity, Phys. Rev. D 8 , 3308 (1973).\n- [42] Y. Hu, P.-P. Wang, Y.-J. Tan, and C.-G. Shao, Testing the polarization of gravitational-wave background with the lisa-tianqin network, The Astrophysical Journal 961 , 116 (2024).\n- [43] S. Capozziello, C. Corda, and M. F. De Laurentis, Massive gravitational waves from f(R) theories of gravity: Potential detection with LISA, Phys. Lett. B 669 , 255 (2008), arXiv:0812.2272 [astro-ph].\n- [44] J. Naf and P. Jetzer, On Gravitational Radiation in Quadratic f ( R ) Gravity, Phys. Rev. D 84 , 024027 (2011), arXiv:1104.2200 [gr-qc].\n- [45] O. Bertolami and A. Martins, On the dynamics of perfect fluids in non-minimally coupled gravity, Phys. Rev. D 85 , 024012 (2012), arXiv:1110.2379 [gr-qc].\n- [46] L. Amendola et al. , Cosmology and fundamental physics with the Euclid satellite, Living Rev. Rel. 21 , 2 (2018), arXiv:1606.00180 [astro-ph.CO].\n- [47] T. Callister, A. S. Biscoveanu, N. Christensen, M. Isi, A. Matas, O. Minazzoli, T. Regimbau, M. Sakellariadou, J. Tasson, and E. Thrane, Polarization-based tests of gravity with the stochastic gravitational-wave background, Phys. Rev. X 7 , 041058 (2017).\n- [48] C. M. Will, The Confrontation between General Relativity and Experiment, Living Rev. Rel. 17 , 4 (2014), arXiv:1403.7377 [gr-qc].\n- [49] A. Nishizawa, A. Taruya, K. Hayama, S. Kawamura, and M.-a. Sakagami, Probing nontensorial polarizations of stochastic gravitational-wave backgrounds with groundbased laser interferometers, Phys. Rev. D 79 , 082002 (2009).\n- [50] E. Berti, V. Cardoso, and C. M. Will, On gravitationalwave spectroscopy of massive black holes with the space\n- interferometer LISA, Phys. Rev. D 73 , 064030 (2006), arXiv:gr-qc/0512160.\n- [51] J. Crowder and N. J. Cornish, Beyond LISA: Exploring future gravitational wave missions, Phys. Rev. D 72 , 083005 (2005), arXiv:gr-qc/0506015.\n- [52] L. Philippoz and P. Jetzer, Detecting additional polarization modes with LISA, J. Phys. Conf. Ser. 840 , 012057 (2017).\n- [53] A. Nishizawa, A. Taruya, and S. Kawamura, Cosmological test of gravity with polarizations of stochastic gravitational waves around 0.1-1 hz, Phys. Rev. D 81 , 104043 (2010).\n- [54] A. Mantziris and O. Bertolami, Gravitational waves from a curvature-induced phase transition of a Higgsportal dark matter sector, arXiv:2407.18845 [astroph.CO] (2024).\n- [55] A. Nishizawa, S. Kawamura, T. Akutsu, K. Arai, K. Yamamoto, D. Tatsumi, E. Nishida, M.-a. Sakagami, T. Chiba, R. Takahashi, and N. Sugiyama, Laserinterferometric detectors for gravitational wave backgrounds at 100 mhz: Detector design and sensitivity, Phys. Rev. D 77 , 022002 (2008).\n- [56] N. Aggarwal et al. , Challenges and opportunities of gravitational-wave searches at MHz to GHz frequencies, Living Rev. Rel. 24 , 4 (2021), arXiv:2011.12414 [gr-qc].\n- [57] P. Navarro, B. Gimeno, J. Monz'o-Cabrera, A. D'ıazMorcillo, and D. Blas, Study of a cubic cavity resonator for gravitational waves detection in the microwave frequency range, Phys. Rev. D 109 , 104048 (2024).\n- [58] C. Gatti, L. Visinelli, and M. Zantedeschi, Cavity detection of gravitational waves: Where do we stand?, Phys. Rev. D 110 , 023018 (2024).\n- [59] F.-Y. Li, H. Wen, Z.-Y. Fang, D. Li, and T.-J. Zhang, Electromagnetic response to high-frequency gravitational waves having additional polarization states: distinguishing and probing tensor-mode, vector-mode and scalar-mode gravitons, Eur. Phys. J. C 80 , 879 (2020), arXiv:1712.00766 [gr-qc].\n- [60] M. Kramer, D. C. Backer, J. M. Cordes, T. J. W. Lazio, B. W. Stappers, and S. Johnston, Strong-field tests of gravity using pulsars and black holes, New Astron. Rev. 48 , 993 (2004), arXiv:astro-ph/0409379.\n- [61] M. E. d. S. Alves and M. Tinto, Pulsar timing sensitivities to gravitational waves from relativistic metric theories of gravity, Phys. Rev. D 83 , 123529 (2011).\n- [62] W. Qin, K. K. Boddy, M. Kamionkowski, and L. Dai, Pulsar-timing arrays, astrometry, and gravitational  \n- waves, Phys. Rev. D 99 , 063002 (2019).\n- [63] Y.-M. Wu, Z.-C. Chen, and Q.-G. Huang, Constraining the polarization of gravitational waves with the parkes pulsar timing array second data release, The Astrophysical Journal 925 , 37 (2022).\n- [64] J. R. Gair, J. D. Romano, and S. R. Taylor, Mapping gravitational-wave backgrounds of arbitrary polarisation using pulsar timing arrays, Phys. Rev. D 92 , 102003 (2015).\n- [65] K. J. Lee, F. A. Jenet, and R. H. Price, Pulsar timing as a probe of non-einsteinian polarizations of gravitational waves, The Astrophysical Journal 685 , 1304 (2008).\n- [66] V. Salvatelli, F. Piazza, and C. Marinoni, Constraints on modified gravity from Planck 2015: when the health of your theory makes the difference, JCAP 09 , 027, arXiv:1602.08283 [astro-ph.CO].\n- [67] S. F. Daniel and E. V. Linder, Confronting General Relativity with Further Cosmological Data, Phys. Rev. D 82 , 103523 (2010), arXiv:1008.0397 [astro-ph.CO].\n- [68] C. Garcia-Quintero, M. Ishak, and O. Ning, Current constraints on deviations from General Relativity using binning in redshift and scale, JCAP 12 , 018, arXiv:2010.12519 [astro-ph.CO].\n- [69] E.-M. Mueller, W. Percival, E. Linder, S. Alam, G.-B. Zhao, A. G. S'anchez, F. Beutler, and J. Brinkmann, The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: constraining modified gravity, Mon. Not. Roy. Astron. Soc. 475 , 2122 (2018), arXiv:1612.00812 [astro-ph.CO].\n- [70] T. M. C. Abbott et al. (DES), Dark Energy Survey Year 3 results: Constraints on extensions to ΛCDM with weak lensing and galaxy clustering, Phys. Rev. D 107 , 083504 (2023), arXiv:2207.05766 [astro-ph.CO].\n- [71] N. Aghanim et al. (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 , A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO].\n- [72] O. Bertolami and J. P'aramos, On the non-trivial gravitational coupling to matter, Class. Quant. Grav. 25 , 245017 (2008), arXiv:0805.1241 [gr-qc].\n- [73] R. van de Weygaert and E. Platen, Cosmic Voids: structure, dynamics and galaxies, Int. J. Mod. Phys. Conf. Ser. 01 , 41 (2011), arXiv:0912.2997 [astro-ph.CO].\n- [74] K. Abazajian et al. (CMB-S4), CMB-S4: Forecasting Constraints on Primordial Gravitational Waves, Astrophys. J. 926 , 54 (2022), arXiv:2008.12619 [astro-ph.CO]."}