Resonant Decay of Gravitational Waves into Dark Energy

TL;DR

This paper investigates how gravitational waves decay into dark energy fluctuations within scalar-tensor theories, revealing potential observable effects on gravitational wave signals detectable by LIGO/Virgo and LISA.

Contribution

It provides an analytical study of gravitational wave decay into dark energy fluctuations, highlighting the impact of specific operators and the potential for observational detection.

Findings

Resonance effects can significantly modify gravitational wave signals.

Observable modifications are predicted within specific parameter ranges for current and future detectors.

Self-interactions of dark energy fluctuations influence the instability and decay process.

Abstract

We study the decay of gravitational waves into dark energy fluctuations $\pi$, taking into account the large occupation numbers. We describe dark energy using the effective field theory approach, in the context of generalized scalar-tensor theories. When the $m_3^3$ (cubic Horndeski) and $\tilde m_4^2$ (beyond Horndeski) operators are present, the gravitational wave acts as a classical background for $\pi$ and modifies its dynamics. In particular, $\pi$ fluctuations are described by a Mathieu equation and feature instability bands that grow exponentially. Focusing on the regime of small gravitational-wave amplitude, corresponding to narrow resonance, we calculate analytically the produced $\pi$, its energy and the change of the gravitational-wave signal. The resonance is affected by $\pi$ self-interactions in a way that we cannot describe analytically. This effect is very relevant for the operator $m_3^3$ and it limits the instability. In the case of the $\tilde m_4^2$ operator self-interactions can be neglected, at least in some regimes. The modification of the gravitational-wave signal is observable for $3 \times 10^{-20} \lesssim \alpha_{\rm H} \lesssim 10^{-17}$ with a LIGO/Virgo-like interferometer and for $10^{-16} \lesssim \alpha_{\rm H} \lesssim 10^{-10}$ with a LISA-like one.

Figures

Equations

$$S=S_0+S_{m_3}+S_{{\tilde{m}}_4},$$

$$S_0$$

$$S_{m_3}$$

$$S_{{\tilde{m}}_4}$$

$$g^{00}$$

$$\delta K$$

$$\text{d}{s}^2=-(1+2\Phi )\text{d}{t}^2+a^2(t)(1-2\Psi )(e^{\gamma }{)}_{ij}\text{d}{x}^i\text{d}{x}^j,$$

$$2M_{\mathrm{Pl}}^2{\nabla }^2\Psi +m_3^3{\nabla }^2\pi =0,M_{\mathrm{Pl}}^2{\nabla }^2(\Phi -\Psi )=0,$$

$$\Phi =\Psi =-\frac{m_3^3}{2M_{\mathrm{Pl}}^2}\pi .$$

$$\alpha \equiv \frac{4M_{\mathrm{Pl}}^2(c+2m_2^4)+3m_3^6}{2M_{\mathrm{Pl}}^4{H}^2},$$

$${\pi }_c\equiv \sqrt{\alpha }{M}_{\mathrm{Pl}}H\pi ,{\gamma }_{ij}^c\equiv \frac{M_{\mathrm{Pl}}}{\sqrt{2}}{\gamma }_{ij}.$$

$$-\frac{1}{2}\sqrt{-g}{m}_3^3\delta K\delta g^{00}\supset \frac{1}{{\Lambda }^2}{\dot{\gamma }}_{ij}^c{\partial }_i{\pi }_c{\partial }_j{\pi }_c,$$

$${\Lambda }^2\equiv \frac{4M_{\mathrm{Pl}}^2(c+2m_2^4)+3m_3^6}{\sqrt{2}{m}_3^3{M}_{\mathrm{Pl}}}=-\frac{\alpha }{\sqrt{2}{\alpha }_{\mathrm{B}}}\frac{H}{H_0}{\Lambda }_2^2,$$

$${\alpha }_{\mathrm{B}}-\frac{{\alpha }_{\mathrm{M}}}{2}(1-{\beta }_1)+{\beta }_1-{\dot{\beta }}_1/H.$$

$${\alpha }_{\mathrm{B}}\equiv -\frac{m_3^3}{2M_{\mathrm{Pl}}^{2}H}.$$

$$S_0+S_{m_3}=\int {\text{d}}^{4}x\left[\frac{1}{4}({\dot{\gamma }}_{ij}{)}^2-\frac{1}{4}({\partial }_k{\gamma }_{ij}{)}^2+\frac{1}{2}{\dot{\pi }}^2-\frac{c_s^2}{2}({\partial }_i\pi {)}^2+\frac{1}{{\Lambda }^2}{\dot{\gamma }}_{ij}{\partial }_i\pi {\partial }_j\pi \right],$$

$$c_s^2=\frac{4M_{\mathrm{Pl}}^{2}c-m_3^3(m_3^3-2M_{\mathrm{Pl}}^{2}H)}{4M_{\mathrm{Pl}}^2(c+2m_2^4)+3m_3^6}=\frac{2}{\alpha }\left(\frac{c}{M_{\mathrm{Pl}}^2{H}^2}-{\alpha }_{\mathrm{B}}-{\alpha }_{\mathrm{B}}^2\right).$$

$${\Gamma }_{\gamma \to \pi \pi }=\frac{p^5(1-c_s^2{)}^2}{480\pi c_s^7{\Lambda }^4},$$

$$\ddot{\pi }-c_s^2{\nabla }^2\pi +\frac{2}{{\Lambda }^2}{\dot{\gamma }}_{ij}{\partial }_i{\partial }_j\pi =0.$$

$${\gamma }_{ij}=M_{\mathrm{Pl}}{h}^{+}{ϵ}_{ij}^{+},h^{+}(t,z)\equiv h_0^{+}\sin (\omega (t-z)),$$

$$\ddot{\pi }-c_s^2{\nabla }^2\pi +c_s^2\beta \cos [\omega (t-z)]({\partial }_x^2-{\partial }_y^2)\pi =0,$$

$$\beta \equiv \frac{2\omega M_{\mathrm{Pl}}{h}_0^{+}}{c_s^2|{\Lambda }^2|}=\frac{2\sqrt{2}|{\alpha }_{\mathrm{B}}|}{\alpha c_s^2}\frac{\omega }{H}{h}_0^{+},\text{for}{m}_3^3\ne 0,{\tilde{m}}_4^2=0({\alpha }_{\mathrm{B}}\ne 0,{\alpha }_{\mathrm{H}}=0).$$

$$S_0+S_{{\tilde{m}}_4}=\int {\text{d}}^{4}x\left[\frac{1}{4}({\dot{\gamma }}_{ij}{)}^2-\frac{1}{4}({\partial }_k{\gamma }_{ij}{)}^2+\frac{1}{2}{\dot{\pi }}^2-\frac{c_s^2}{2}({\partial }_i\pi {)}^2+\frac{1}{{\Lambda }_{\star }^3}{\ddot{\gamma }}_{ij}{\partial }_i\pi {\partial }_j\pi \right],$$

$$c_s^2=\frac{2}{\alpha }\left[(1+{\alpha }_{\mathrm{H}}{)}^2\frac{c}{M_{\mathrm{Pl}}^2{H}^2}+{\alpha }_{\mathrm{H}}+{\alpha }_{\mathrm{H}}(1+{\alpha }_{\mathrm{H}})\frac{\dot{H}}{H^2}\right],$$

$${\Lambda }_{\star }^3\equiv \frac{\sqrt{2}{M}_{\mathrm{Pl}}^3(c+2m_2^4)}{{\widetilde{m_4}}^2(M_{\mathrm{Pl}}^2+2{\tilde{m}}_4^2)}\simeq \sqrt{2}\frac{\alpha }{{\alpha }_{\text{H}}}{\left(\frac{H}{H_0}\right)}^2{\Lambda }_3^3,$$

$$\ddot{\pi }-c_s^2{\nabla }^2\pi +\frac{2}{{\Lambda }_{\star }^3}{\ddot{\gamma }}_{ij}{\partial }_i{\partial }_j\pi =0.$$

$$\ddot{\pi }-c_s^2{\nabla }^2\pi -\frac{2{\omega }^2{M}_{\mathrm{Pl}}{h}_0^{+}}{{\Lambda }_{\star }^3}\sin [\omega (t-z)]({\partial }_x^2-{\partial }_y^2)\pi =0.$$

$$\beta \equiv \frac{2{\omega }^2{M}_{\mathrm{Pl}}{h}_0^{+}}{c_s^2|{\Lambda }_{\star }^3|}=\frac{\sqrt{2}|{\alpha }_{\mathrm{H}}|}{\alpha c_s^2}{\left(\frac{\omega }{H}\right)}^2{h}_0^{+},\text{for}{m}_3^3=0,{\tilde{m}}_4^2\ne 0({\alpha }_{\mathrm{B}}=0,{\alpha }_{\mathrm{H}}\ne 0).$$

$$u\equiv t-z,v\equiv t+z,$$

$${\gamma }_{ij}(u)=M_{\mathrm{Pl}}{h}_0^{+}\sin \left(\omega u\right){ϵ}_{ij}^{+},$$

$$\ddot{\pi }-c_s^2{\nabla }^2\pi +c_s^2\beta \cos \left(\omega u\right)({\partial }_x^2-{\partial }_y^2)\pi =0.$$

$$s\equiv -t+c_s^{-2}z,$$

Full text

1. 1 Introduction

2. 2 Graviton-scalar-scalar vertices

3. 2.1 $m_{3}^{3}$-operator

4. 2.2 $\tilde{m}_{4}^{2}$-operator

5. 3 Narrow resonance

6. 3.1 Parametric resonance

7. 3.2 Energy density of $\pi$

8. 3.3 Modification of the gravitational waveform

9. 3.4 Generic polarization

10. 3.5 Conservation of energy

11. 4 Nonlinearities

12. 4.1 $m_{3}^{3}$-operator

13. 4.2 $\tilde{m}_{4}^{2}$-operator

14. 5 Observational signatures for $\tilde{m}_{4}^{2}$

15. 5.1 Fundamental frequency

16. 5.2 Higher harmonics and precursors

17. 6 Conclusions and future directions

18. A Interactions in Spatially-Flat Gauge

19. B Parametric resonance as Bose enhancement

20. C Details on the conservation of energy

21. C.1 $m_{3}^{3}$-operator

22. C.2 $\tilde{m}_{4}^{2}$-operator

1 Introduction

In a recent article [1], some of us pointed out that in modified-gravity theories gravitational waves (GWs) can decay into dark energy fluctuations, as a consequence of the spontaneous breaking of Lorentz invariance. We focused on dark energy and modified gravity based on a scalar degree of freedom, whose time-dependent background induces a preferred time-foliation of the FRW metric. The natural way to describe this setting is the effective field theory of dark energy (EFT of DE) [2, 3, 4, 5, 6], which we review below. Moreover, since the relative speed between GWs and light is now constrained with $10^{-15}$ accuracy [7], we considered only models where gravitons travel at the same speed as light [8, 9, 10, 11]. (One possible caveat about these constraints is that they do not apply if the theory breaks down at scales parametrically lower than the frequency of the observed GWs [12]. In this case the propagation of GWs, as well as many accurate tests of gravity, can only be described if one knows the UV completion of the EFT of DE. )

In particular, we showed that some of the operators of the EFT of DE display a cubic $\gamma\pi\pi$ interaction, where $\gamma$ and $\pi$ respectively denote the graviton and the scalar-field fluctuation, that can mediate the decay. Of course, this can happen only if both energy and momentum can remain conserved during the process, which is when scalar fluctuations propagate subluminally. The same vertex is also responsible for an anomalous GW dispersion, for speeds of scalar propagation different from that of light. Depending on the energy scale suppressing this interaction, these two effects can be important at frequencies observed by LIGO/Virgo and can constrain these theories.

The bound derived in [1] is based on a perturbative calculation, in which individual gravitons are assumed to decay independently of each other. But a classical GW is a collection of many particles with very large occupation number and particle production must be treated as a collective process in which many gravitons decay simultaneously. The classical GW acts as a background for the propagation of scalar fluctuations. This paper studies this process in the limit where the GW background acts as a small periodic perturbation (narrow resonance). Larger amplitudes of the GW can induce tachyon or ghost instabilities: we will study this possibility in a forthcoming publication [13].

To review the EFT approach, we define the necessary geometrical quantities. We assume the flat FRW metric to be given by $\text{d}s^{2}=-\text{d}t^{2}+a^{2}(t)\text{d}\boldsymbol{x}^{2}$, so that $\delta g^{00}\equiv 1+g^{00}$ is the perturbation of $g^{00}$ around the background solution. Moreover, we define by $\delta K^{\mu}_{\ \nu}\equiv K^{\mu}_{\ \nu}-H\delta^{\mu}_{\ \nu}$, where $H\equiv\dot{a}/a$ is the Hubble rate, the perturbation of the extrinsic curvature of the equal-time hypersurfaces and by $\delta K$ its trace. Finally, we will need the 3d Ricci scalar of these hypersurfaces, ${}^{(3)}\!R$.

As mentioned above, we focus on theories where gravitons travel luminally and, for later convenience, we split the EFT of DE action in the sum of three actions [8],

[TABLE]

where

[TABLE]

The first action, $S_{0}$, contains the Einstein-Hilbert term and the minimal scalar field Lagrangian, which describe the dynamics of the background. Notice that we have removed any time-dependence in front of the Einstein-Hilbert term by a conformal transformation, which leaves the graviton speed unaffected. Indeed, the first two functions of time, $c(t)$ and $\lambda(t)$,111We use the notation $\lambda(t)$, and not $\Lambda(t)$ as usual in the literature, to avoid confusion with the energy scale $\Lambda$ suppressing the cubic vertices studied in this article. can be determined by the background evolution in terms of the Hubble expansion and matter quantities (see e.g. Sec. 5 of [5] for their explicit expressions in this frame). The operator $m_{2}^{4}$ does not change the background but affects the speed of propagation of scalar fluctuations, $c_{s}^{2}$. Its typical value is ${\cal O}(M_{\rm Pl}^{2}H^{2})$, see e.g. [5] for details. This action was studied in details in [4] and in the covariant language it describes quintessence [14] or, more generally, a dark energy with scalar field Lagrangian $P(\phi,X)$ [15], where $X\equiv g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi$.

The operator in the second action, $S_{m_{3}}$, introduces a kinetic mixing between the scalar field and gravity [2, 4] (sometimes called kinetic gravity braiding [16, 17]). In the covariant language it corresponds to the cubic Horndeski Lagrangian [18, 19], of the form $Q(\phi,X)\Box\phi$ with $\Box\equiv g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$. In the regime that leads to sizeable modifications of gravity (i.e. for $m_{3}^{3}\sim M_{\rm Pl}^{2}H_{0}$), the operator contained in $S_{m_{3}}$ displays a $\gamma\pi\pi$ interaction suppressed by an energy scale of order $\Lambda_{2}\equiv(M_{\rm Pl}H_{0})^{1/2}\sim 10^{-3}$eV. This energy scale is much greater than the typical LIGO/Virgo frequency. For this reason, in [1] the parameter $m_{3}^{3}$ remains unconstrained by the graviton decay computed in perturbation theory. Finally, the operator in the third action, $S_{\tilde{m}_{4}}$ whose typical size is $M_{\rm Pl}H_{0}$ introduces a kinetic mixing between the scalar field and matter [20] (or kinetic matter mixing [21]). In the covariant language it corresponds to the beyond-Horndeski Lagrangian of the Gleyzes-Langlois-Piazza-Vernizzi type [20, 22]. This operator displays a $\gamma\pi\pi$ interaction and was constrained in [1] with the perturbative decay because the vertex is suppressed by an energy scale close to LIGO/Virgo frequencies.

In Sec. 2, after expanding the action $S_{0}+S_{m_{3}}$ in perturbations in the Newtonian gauge (the same calculation is repeated in the spatially-flat gauge in App. A), we study the effect of a classical GW background on the $\pi$ dynamics, for the operator $m_{3}^{3}$ (Sec. 2.1) and $\tilde{m}_{4}^{2}$ (Sec. 2.2). The regime of small GWs can be studied analytically and leads to the so-called narrow parametric resonance, which is the subject of Sec. 3. There we compute the energy density of $\pi$ produced by the parametric instability due to the oscillating GWs (in Sec. 3.2) and we re-interpret the $\pi$ production in the narrow-resonance regime as an effect of Bose enhancement of the perturbative decay in App. B. The back-reaction on the GW signal is computed in Sec. 3.3 for a linearly polarized wave, while the case of elliptical polarization is discussed in Sec. 3.4. In Sec. 3.5 we check that energy is conserved in this process, as expected (the details of the calculations are given in App. C).

The treatment in Sec. 3 neglects scalar-field nonlinearities, which are studied in Sec. 4. The operator $m_{3}^{3}$ contains cubic self-interactions suppressed by the scale $\Lambda_{3}\equiv(M_{\rm Pl}H_{0}^{2})^{1/3}\sim 10^{-13}$eV, which is much smaller than the one appearing in the vertex $\gamma\pi\pi$. Thus, these become relevant and probably halt the parametric resonance well before the GWs are affected by the back-reaction (Sec. 4.1). This makes the results of Sec. 3 applied to this operator inconclusive. The situation is different for the operator $\tilde{m}_{4}^{2}$: in this case the scale that suppresses non-linearities is the same that appears in the coupling $\gamma\pi\pi$. The leading non-linearities are quartic in the regime of interest and are suppressed with respect to a naive estimate due to the particular structure of Galileon interactions. At least in some region of parameters non-linearities do not halt the parametric instability due to the oscillating GWs. In Sec. 5.1 we therefore study in which range of parameters one expects a modification of the GW signal. Moreover, in Sec. 5.2 we discuss precursors, higher harmonics induced in the GW signal by the produced $\pi$, that enter the observational band earlier than the main signal. We conclude discussing the main results of the article and possible future directions in Sec. 6.

2 Graviton-scalar-scalar vertices

Let us derive the interaction $\gamma\pi\pi$ from the action (1.1), using the Newtonian gauge. For the time being, we neglect the self-interactions of the $\pi$ field; they will be discussed later, in Sec. 4. We initially focus on the operator $m_{3}^{3}$; as a check, in App. A we perform the same calculation in spatially-flat gauge.

2.1 $m_{3}^{3}$-operator

Let us consider the action $S_{0}+S_{m_{3}}$. One can restore the $\pi$ dependence in a generic gauge with the Stuekelberg procedure [3, 6] $t\rightarrow t+\pi(t,\boldsymbol{x})$. Focusing on the terms relevant for our calculations we have [23]

[TABLE]

In Newtonian gauge, the line element reads

[TABLE]

with $\gamma$ transverse, $\partial_{i}\gamma_{ij}=0$, and traceless, $\gamma_{ii}=0$.

Varying the above action with respect to $\Phi$ and $\Psi$ and focussing on the sub-Hubble limit by keeping only the leading terms in spatial derivatives, one obtains

[TABLE]

where $\nabla^{2}\equiv\partial_{x}^{2}+\partial_{y}^{2}+\partial_{z}^{2}$. From now on we will always consider the Minkowski limit, i.e. that time and spatial derivatives are much larger than Hubble. These equations can be solved in terms of $\pi$,

[TABLE]

Using these relations, one can find the kinetic term of the $\pi$ field.

To write this, it is convenient to define the dimensionless quantity [24, 25]

[TABLE]

which sets the normalization of the scalar fluctuations (and must be positive to avoid ghost instabilities). The canonically normalised scalar and tensor perturbations are then given by

[TABLE]

In terms of these, the interaction term, after integrating by parts, reads

[TABLE]

where

[TABLE]

and in the right-hand side we have defined the dimensionless quantity222To write the action (1.1), we used a conformal transformation (possibly dependent on $X=(\partial_{\mu}\phi)^{2}$) to set to constant the effective Planck mass and to zero higher-order operators of DHOST theories [26, 27]. Using the notation of [28] and the transformation formulae contained therein, one can check that in a general frame the only relevant parameter is

$\alpha_{\rm B}-\frac{\alpha_{\rm M}}{2}(1-\beta_{1})+\beta_{1}-\dot{\beta}_{1}/H\;.$

(2.10)

[TABLE]

We drop the symbol of canonical normalisation: $\gamma$ and $\pi$ will indicate for the rest of the paper the canonically normalised fields. The total Lagrangian is then

[TABLE]

and

[TABLE]

The perturbative decay rate of the graviton can be computed following exactly the same procedure followed for $\tilde{m}_{4}^{2}$ in [1]. It gives

[TABLE]

where $p$ is the momentum of the decaying graviton. One can check that this is negligible for frequencies relevant for GW observations, since by eq. (2.9) $\Lambda$ is of order $\Lambda_{2}$.

The equation of motion of $\pi$ from the Lagrangian (2.12) reads

[TABLE]

Let us use the classical background solution of the GW travelling in the $\hat{z}$ direction with a linear polarization. Without loss of generality we take the $+$ polarization (the $\times$ one can be obtained by a $45^{\rm o}$ rotation of the axes)

[TABLE]

where $h^{+}$ is the dimensionless strain of the GW and the polarisation tensor is $\epsilon_{ij}^{+}={\rm diag}(1,-1,0)$. In this paper we will always be away from the source generating the GW, i.e. in the weak field regime $h^{+}\ll 1$. Substituting the solution (2.16) into eq. (2.15), one gets (modulo an irrelevant phase)

[TABLE]

where the parameter $\beta$ is defined as

[TABLE]

In the following, to evaluate the right-hand side of the above definition we will use $\omega=\Lambda_{3}$ and $H=H_{0}$, in which case $\omega/H\sim 10^{20}$. Moreover, note that $\alpha$ and $c_{s}^{2}$ in the second equality appear in the combination $\alpha c_{s}^{2}$. Therefore, thanks to eqs. (2.13), $\beta$ defined above depends only on $\alpha_{\rm B}$ (or $m_{3}$) and we can tune this parameter to make $\beta$ small.

2.2 $\tilde{m}_{4}^{2}$-operator

We now consider the operator $\tilde{m}_{4}^{2}$. To avoid large $\pi$ nonlinearities, discussed in Sec. 4, we focus on the case $m_{3}=0$ and consider the action $S_{0}+S_{\tilde{m}_{4}}$. The operator $\tilde{m}_{4}^{2}$ has been studied in [1] and details on the calculations can be found there.

In this case, the action for the canonically normalised fields reads

[TABLE]

where the sound speed of scalar fluctuations is now

[TABLE]

and we have defined $\alpha_{\rm H}\equiv 2\tilde{m}_{4}^{2}/M_{\rm Pl}^{2}$ (we assumed $\alpha_{\rm H}\ll 1$ since this will be the regime of interest). Thus, this operator does not affect the speed of propagation of GWs [8], but, as shown in [1], it contains an interaction $\gamma\pi\pi$ suppressed by the scale

[TABLE]

where in the last equation we assumed a small $\alpha_{\rm H}$. This should be compared with the cubic coupling discussed above, eq. (2.12).

The evolution equation for $\pi$ then reads

[TABLE]

Substituting the solution (2.16) into this equation gives

[TABLE]

One sees that the evolution equation for $\pi$ for the $\tilde{m}_{4}^{2}$ operator is very similar to the case of the $m_{3}^{3}$ operator, with the replacement $\Lambda^{2}\to\Lambda_{\star}^{3}\omega^{-1}$. We can thus apply eq. (2.17), with $\beta$ now defined as

[TABLE]

Analogously to the $m_{3}$ case, because of eqs. (2.20), $\beta$ defined above depends only on $\alpha_{\rm H}$ (or $\tilde{m}_{4}$) and also in this case we can tune this parameter to make $\beta$ small.

3 Narrow resonance

Equation (2.17) describes a harmonic oscillator with periodic time-dependent frequency, which can lead to parametric resonance. As explained in the introduction, in this article we are going to focus on the narrow-resonance regime, which corresponds to $\beta\ll 1$. In this case, the solution of the equation of motion of $\pi$ can be treated analytically and features an exponential growth of scalar fluctuations.

3.1 Parametric resonance

The GW is emitted at $t=0$ in the $z$ direction and is detected at some later time, see Fig. 1. Using the light-cone coordinates,

[TABLE]

its solution for $t>0$ can be written as

[TABLE]

where $h_{0}^{+}$ can be taken as constant, since it varies slowly compared to the GW frequency.

Eq. (2.17) takes the form

[TABLE]

For $c_{s}<1$, which we assume in the following, $u$ is a time-like variable for the $\pi$ metric: hypersurfaces of constant $u$ are space-like to the $\pi$-cone. It is then convenient to define the variable $s$,

[TABLE]

which is orthogonal to $u$ and is thus space-like with the $\pi$ metric, see Fig. 1, and use the coordinates $\tilde{\boldsymbol{x}}=(x,y,s)$ to describe the spatial foliations. Since the $\pi$-lightcone is narrower than the one of GWs, the solution for the scalar will only be sensitive to a finite region of the GW background. This is the reason why one can approximate the GW as a plane wave with constant amplitude and disregard the process of emission of the GW from the astrophysical source. In particular, for this to be a good approximation we need to require that the past light-cone of $\pi$ overlaps only with a region of the GW with constant amplitude. From Fig. 1 one sees that this is a very weak requirement. It is enough that $1-c_{s}$ is larger than the ratio between the duration of the GW signal (of order seconds) and the scale of variation of the amplitude (of order Mpc). Plugging the numbers one gets $1-c_{s}\gtrsim 10^{-14}$.

Since for $\pi$ there is translational invariance in $\tilde{\boldsymbol{x}}$, it is useful to decompose $\pi$ in Fourier modes as

[TABLE]

where $\tilde{\boldsymbol{p}}=(p_{x},p_{y},p_{s})$ is conjugate to $\tilde{\boldsymbol{x}}$. In the absence of the interaction with the gravitational wave, one can relate $p_{s}$ to the momentum written in the original coordinates,

[TABLE]

In the following we are going to use this change of variable also when $\beta\neq 0$, although plane waves in the original coordinates $(t,x,y,z)$ are not solution of eq. (3.3).

Then we can quantize $\pi$ straightforwardly. More specifically, we decompose $\pi_{\tilde{\boldsymbol{p}}}$ as

[TABLE]

where

[TABLE]

and $\hat{a}_{\tilde{\boldsymbol{p}}}$ and $\hat{a}_{-\tilde{\boldsymbol{p}}}^{\dagger}$ are the usual creation and annihilation operators satisfying the commutation relations, $[\hat{a}_{\tilde{\boldsymbol{p}}},\hat{a}^{\dagger}_{\tilde{\boldsymbol{p}}^{\prime}}]=(2\pi)^{3}\delta^{(3)}(\tilde{\boldsymbol{p}}-\tilde{\boldsymbol{p}}^{\prime})$. The normalization is chosen for convenience. Indeed, for $\beta=0$ the evolution equation for $\pi$ satisfies a free wave equation and each Fourier mode can be described as an independent quantum harmonic oscillator. We assume that in this case $\pi$ is in the standard Minkowski vacuum, given by333It is straightforward to verify that eq. (3.10) is equivalent to the standard Minkowski vacuum, i.e.

$\pi(x)=\int\frac{\text{d}^{3}{\boldsymbol{p}}}{(2\pi)^{3}}\frac{1}{\sqrt{2c_{s}|{\boldsymbol{p}}|}}\left(e^{-ip\cdot x}\hat{a}_{\boldsymbol{p}}+e^{ip\cdot x}\hat{a}^{\dagger}_{-\boldsymbol{p}}\right)\;,$

(3.9) upon use of $\text{d}^{3}\tilde{\boldsymbol{p}}/\text{d}^{3}\boldsymbol{p}={c_{s}p_{u}}/{|\boldsymbol{p}|}$ and, consequently, of $\hat{a}_{\tilde{\boldsymbol{p}}}=\left[{|{\boldsymbol{p}}|}/({c_{s}p_{u}})\right]^{1/2}\hat{a}_{{\boldsymbol{p}}}$.

[TABLE]

To study the parametric resonance, will now show that eq. (3.3) can be written as a Mathieu equation [29]. First, in terms of the new coordinates, eq. (3.3) becomes

[TABLE]

For convenience we can also define the dimensionless time variable $\tau$,

[TABLE]

For each Fourier mode, $f$ satisfies

[TABLE]

with

[TABLE]

To write $A$ and $q$ we have decomposed the vector $\boldsymbol{p}$ in polar coordinates, $\boldsymbol{p}=|\boldsymbol{p}|(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta)$, and we have defined $\Omega\equiv p_{z}/|\boldsymbol{p}|=\cos\theta$.

The general solution of the Mathieu equation is of the form $e^{\pm\mu\tau}P(\tau)$, where $P(\tau+\pi)=P(\tau)$ [29]. If the characteristic exponent $\mu$ has a real part, the solution of the Mathieu equation is unstable for generic initial conditions. Since $\mu$ is a function of $A$ and $q$, the instability region can be represented on the $(q,A)$-plane, see Fig. 2(a). The unstable regions are also shown in Fig. 2(b) on the plane $(\Omega,{c_{s}^{2}\,\boldsymbol{p}^{2}}/{\omega^{2}})$, for specific values of the other parameters (in the example in the figure we take $\beta=0.8$ and $c_{s}=1/2$). Notice that the maximal exponential growth is reached when the ratio $q/A$ has its maximum at $\Omega=c_{s}$ and $\varphi=0$.

3.2 Energy density of $\pi$

In this subsection we compute the energy density of $\pi$ produced by the coupling with the GW. This is given by (see eq. (3.52) and explanation below)

[TABLE]

Decomposing in Fourier modes and in the $\hat{a}$ and $\hat{a}^{\dagger}$ coefficients using respectively eqs. (3.5) and (3.7), the energy density can be rewritten as

[TABLE]

where we have simplified the expression on the right-hand side using that the Wronskian is time-independent, $\mathscr{W}[f_{\tilde{p}}(u),f^{\star}_{\tilde{p}}(u)]=-2ip_{u}$.444In general, we can write $f$ as a linear combination of Mathieu-sine and cosine functions, respectively $\mathcal{S}$ and $\mathcal{C}$ [29]. We fix the boundary conditions of the solution such that $\pi$ is in the vacuum, i.e. the function $f$ satisfies eq. (3.10), at $u=0$. This gives

$f(\tau)=\frac{-i2p_{u}}{\omega}\,\mathcal{S}\left(A,q;\tau\right)+\,\mathcal{C}\left(A,q;\tau\right)\;.$

(3.18)

Using this expression and the properties of these functions, we can check that the Wronskian $\mathscr{W}[f_{\tilde{p}}(u),f^{\star}_{\tilde{p}}(u)]$ is constant and with the above normalization is given by $-2ip_{u}$. As a consequence, the commutation relation in the “interacting” region are satisfied at all times.

One can verify that in the limit $\beta=0$ the above expression reduces to the energy density of the vacuum, i.e. $\rho_{\pi}^{0}=\int\frac{\text{d}^{3}p}{(2\pi)^{3}}\frac{\omega_{p}}{2}$.

We can simplify the right-hand side further by making some approximations. Since we are not interested in following the oscillatory behaviour of $\rho_{\pi}$, we can perform an average in $\tau$ over many periods of oscillation. Since the amplitude of the periodic part of $f_{\tilde{\boldsymbol{p}}}$ is bounded to be less than unity, it is reasonable to take $\braket{}{f_{\tilde{p}}}{{}^{2}}_{T}\simeq\braket{}{\partial_{\tau}f_{\tilde{p}}}{{}^{2}}_{T}\simeq e^{2\mu\tau}/2$ in eq. (3.2). This educated guess will be confirmed in Sec. 3.3. Changing integration variables, from $\tilde{\boldsymbol{p}}$ to $\xi\equiv c_{s}^{2}\boldsymbol{p}^{2}/\omega^{2}$ and the angular variables $\Omega$ and $\varphi$, we find

[TABLE]

We now want to solve the integral on the right-hand side. Since it is dominated by the unstable modes, we restrict the domain of integration to the bands of instability. Actually, we are going to focus on the first instability band for two reasons. The instability rate of the higher bands goes as $\mu_{m}\sim q^{m}/(m!)^{2}\sim\beta^{m}$ [29, 30], which implies that for small $\beta$ the first band is the most unstable. Second, the $\pi$ produced in this band will modify the GW signal with the original angular frequency $\omega$. (With some contributions at higher frequencies that will be discussed in Section 5.2.) We are going to work in the regime $\beta\ll 1$ corresponding to $q\ll 1$, the regime of narrow resonance.

In this situation we can restrict the integral to the first unstable band, which is defined by [29]

[TABLE]

Within this region, the value of the exponent $\mu$ is

[TABLE]

Using the definitions (3.14) and (3.15) in eq. (3.20), the boundary region above can be rewritten in terms of $\xi$,

[TABLE]

which fixes the domain of integration in eq. (3.19).

The integral in (3.19) can be then solved with the saddle-point approximation. In general, an integral of the form

[TABLE]

can be approximated for large $\tau$ by

[TABLE]

where $\boldsymbol{X}_{0}$ is the maximum of the exponent, i.e. the solution of $\partial_{i}h(\boldsymbol{X})=0$ with the Hessian ${\rm H}_{ij}\equiv-\partial_{i}\partial_{j}h(\boldsymbol{X}_{0})$ positive definite. (If there are more than one maximum one should sum over them.)

In our case we have to maximize $\mu$ of eq. (3.21), by requiring that $\partial\mu/\partial\xi$, $\partial\mu/\partial\Omega$ and $\partial\mu/\partial\varphi$ vanish. This happens for

[TABLE]

but we discard the solutions with $n$ odd, because they make $q$ (and the integration region) vanish, see eq. (3.15). There are four relevant saddle points $(\varphi=0,\pi/2,\pi,3\pi/2)$ giving the same value to the integral; without loss of generality we choose $\varphi=0$. The expression of the Hessian matrix at the saddle point is

[TABLE]

Working at leading order in $\beta$ and using the saddle-point approximation with these values in eq. (3.19), the energy density of $\pi$ as a function of $u=t-z$ reads

[TABLE]

We can compare $\rho_{\pi}$ with the energy density of the gravitational wave, which is roughly constant,

[TABLE]

For instance, for $\beta=0.1$ and $c_{s}=1/2$, we get ${\rho_{\pi}}\simeq\rho_{\gamma}$ after $\tau/\pi\simeq 750$ cycles.

The exponential growth studied in this section can also be seen as a consequence of Bose enhancement, see Appendix B.

3.3 Modification of the gravitational waveform

The parametric production of $\pi$ suggests that its back-reaction will modify the background gravitational wave. In this section we estimate this effect remaining in the narrow-resonance limit, i.e. $|q|\ll 1$ ($\beta\ll 1$). Here we focus again on the case of the operator $m_{3}^{3}$.

To compute the back-reaction on $\gamma_{ij}$, we start from the action (2.12). The equation of motion for $\gamma_{ij}$ is

[TABLE]

where, given a direction of propagation of the wave $\boldsymbol{n}$, $\Lambda_{ij,kl}(\boldsymbol{n})$ is the projector into the traceless-transverse gauge, defined by

[TABLE]

We focus again on a wave traveling in the $\hat{\boldsymbol{z}}$ direction. Using light-cone variables the equation above becomes

[TABLE]

We can then split the solution into a homogeneous and a forced one,

[TABLE]

The former reads

[TABLE]

while for the latter, which represents the back-reaction due to $\pi$, we find

[TABLE]

Because it is transverse and traceless, the source $J_{ij}$ can be projected into a plus and cross polarization. We will focus on the plus polarization. In this case, we can proceed as in the previous subsection and rewrite the integral in terms of $\Omega$, $\varphi$ and $\xi$,

[TABLE]

Then, we can evaluate this integral in the narrow-resonance regime with the saddle-point approximation. However, we will now use an approximation for the solution $f_{\tilde{p}}(u)$ that takes into account its oscillatory behaviour. In particular, we will split the integral in two regions: $q>0$ and $q<0$. Since the sign of $q$ is controlled by the angle $\varphi$ through $\cos 2\varphi$, this corresponds to splitting the integration over $\varphi$.

For $q>0$, the solution of the Mathieu equation in the first instability band can be approximated by (see pag. 72 of [29])

[TABLE]

where $\mu>0$ and $\sigma\in(-\frac{\pi}{2},0)$ is a parameter, which depends on $A$ and $q$. It is real inside the instability bands, as shown in the instability chart for the Mathieu equation in Fig. 2(a). More specifically, in the first instability band one has

[TABLE]

The coefficients $c_{+}$ and $c_{-}$ can be fixed by demanding that at $\tau=0$ we recover the vacuum solution.

The case $q<0$ can be obtained by noting that when $q$ changes sign, $\mu$ does it as well while $A$ remains the same. The only way to implement this is to consider the simultaneous change $\sigma\rightarrow\sigma^{\prime}=-\sigma-\frac{\pi}{2}$. By performing these two transformations for $q$ and $\sigma$ on $f(\tau)$ one obtains

[TABLE]

We can now integrate the right-hand side of (3.35) starting from the interval $\varphi\in\left(-\frac{\pi}{4},\frac{\pi}{4}\right)$. Replacing the growing mode solution of eq. (3.36) in the integrand, we obtain

[TABLE]

The $\tau$ dependence in $\sin(\tau-\sigma)$ seems to change the saddle point computed by maximizing $\mu$. However, by rewriting the sine as exponential functions one can check that its effect is simply to add an additional constant to the exponent so that the saddle point remains the same as the one we computed in Sec. 3.2, see eq. (3.25). At the saddle point, eqs. (3.37) and (3.38) give $\tan(2\sigma)\simeq 2/\beta$ while $c_{+}$ and $c_{+}^{\prime}$ can be computed by requiring vacuum initial conditions, eq. (3.10), at $\tau=0$. Working for small $\beta$, this gives

[TABLE]

Then, applying eq. (3.24) to the integral in (3.40), the latter can be solved, giving

[TABLE]

where for $c_{+}$ and $\sigma$ we must use the saddle-point values, eq. (3.41). We can now repeat this exercise for each part of the integral, using the growing mode of either eqs. (3.36) or (3.39) depending on the sign of $q$. Summing them together, we obtain

[TABLE]

Replacing in the solution for $\Delta\gamma_{ij}$, eq. (3.34), the expression of $J_{ij}(u)$ of eq. (3.43), with $\sigma$, $c_{+}$ and $c_{+}^{\prime}$ fixed by the saddle-point approximation from eq. (3.41), the back-reaction on the GWs due to $\pi$ reads

[TABLE]

where we have dropped $J_{ij}(0)$ which is negligible at late time. Thus, the back-reaction grows exponentially in $u$, as expected from the growth of the energy of $\pi$, eq. (3.27), and linearly in $v$. Note that the right-hand side of the above equation diverges for $\beta\rightarrow 0$. This is because it has been obtained using the saddle-point approximation, which assumes that $\beta\omega u$ is large. In Appendix B we check that this result reduces to the perturbative calculation when the occupation number is small enough.

Notice that there is no production of cross polarization. Indeed, the integrand of the source term in eq. (3.31) now contains $2p_{x}p_{y}$ instead of $p_{x}^{2}-p_{y}^{2}$. Since $p_{x}p_{y}\propto\sin(2\varphi)=0$ and the saddle points are such that $\sin(2\varphi)=0$, the cross-polarized waves are not generated by the back-reaction of dark energy fluctuations produced by plus-polarized waves and eq. (3.44) represents the full result.

3.4 Generic polarization

So we far we have been discussing linearly polarized waves. Since the resonant effect is non-linear, one cannot simply superimpose the result for linear polarization in order to get a general polarization. In this subsection we are going to consider a more generic polarization state, that we parametrise as follows

[TABLE]

where $0\leq\alpha<2\pi$ is an angle characterizing the GW polarization. Note that the state of polarization for generic $\alpha$ is elliptical, like the one coming from binary systems [31].

Following the same procedure as in Sec. 3.1, the Mathieu eq. (3.13) becomes

[TABLE]

with

[TABLE]

while $A$ remains the same as before. One needs to shift $\tau\rightarrow\tau+\hat{\theta}/2$ in order to use the same form for the Mathieu solution. Given this change in $q$, one can easily obtain the modified saddle points which are given by

[TABLE]

Thus, the saddle points for $\Omega$ and $\varphi$ are unaffected by the angle $\alpha$. From the saddle point found above one can see that choosing $n$ to be even selects the $+$ polarization, whereas $n$ odd corresponds to $\times$ polarization. The exponent $\mu$ of eq. (3.21) can be evaluated on these saddle points and, at leading order in $\beta$, is given by

[TABLE]

Here one can define $\mu^{+}$ and $\mu^{\times}$ corresponding to $+$ and $\times$ polarizations respectively as

[TABLE]

Several comments are in order at this stage. First, for $0<\alpha<\pi/4$ the $+$ contribution in the initial wave is larger. In this case one has $\mu^{+}>\mu^{\times}$ which means that the $+$ polarization dominates also in the back-reaction for $\Delta\gamma_{ij}$. For $\pi/4<\alpha<\pi/2$ the $\times$ mode dominates. For $\alpha=\pi/4$ both polarization states grow with the same rate, meaning that a circularly polarized wave remains circular. Moreover, by setting $\alpha=0$ ($\alpha=\pi/2$) we recover the results of the previous sections for the case of $+$ ($\times$) polarization.

3.5 Conservation of energy

To check our results and to get a better understanding of the system, it is useful to verify that energy is conserved in the production of the $\pi$ field and the corresponding modification of the GW, $\Delta\gamma$. From the Lagrangian (2.12), one can derive the Noether stress-energy tensor, which is conserved on-shell as a consequence of translational invariance

[TABLE]

(Notice that this $T^{\mu}_{\ \nu}$ will be different from the pseudo stress-energy tensor of GR.) Let us consider the region represented in Fig. 3. Given the symmetries of the system, it is useful to take the left and right boundaries as null, instead of time-like, surfaces.

At first sight, it is somewhat puzzling that while both the original GW and the induced $\pi$ are only displaced as time proceeds (see Figs. 1 and 3), so that their contribution to the total energy does not depend on time, $\Delta\gamma$ grows on later time slices since it is proportional to $v$. What we are going to verify is that the variation of the energy between $\partial\mathcal{V}_{1}$ and $\partial\mathcal{V}_{2}$ due to the change in $\Delta\gamma$ is compensated by a flux of energy across the null boundary $\partial\mathcal{V}_{0}$.555Notice that, while the original GW is described by a wave packet localised in a certain interval of $u$, $\pi$ waves are present at arbitrary large $u$ and thus will always contribute to the flux. (There is no flux across the right null boundary $\partial\mathcal{V}_{0}^{\prime}$ since $\pi$ is in the vacuum for $u=0$.)

The Noether stress-energy tensor of the action (2.12) is given by

[TABLE]

Notice that the total energy of the system is simply the sum of the kinetic energy of $\gamma$ and $\pi$ without a contribution due to interactions. (This is a consequence of the interaction term being linear in $\dot{\gamma}$.) This means that the production of $\pi$ must be compensated by a decrease of the $\gamma$ kinetic energy. Using the splitting in eq. (3.32) and defining $\bar{\rho}_{\gamma}\equiv\frac{1}{4}(\dot{\bar{\gamma}}_{ij})^{2}+\frac{1}{4}(\partial_{k}\bar{\gamma}_{ij})^{2}$, $\rho_{\pi}\equiv\frac{1}{2}\dot{\pi}^{2}+\frac{c_{s}^{2}}{2}(\partial_{i}\pi)^{2}$, the components (3.52) (3.54) can be written up to second order in perturbation. For instance

[TABLE]

Gauss theorem works also when the region has null boundaries (see for instance [32]):

[TABLE]

The only subtlety in the case of null boundaries is that one does not know how to choose the normalization of the null vector $n^{\mu}$ orthogonal to the surface (of course when the boundary is null this vector also lies on the surface). A related ambiguity is that there is no natural volume form on the boundary to perform the integration, since the induced metric on a null surface is degenerate. The two ambiguities compensate each other. If one chooses a 3-form $\tilde{\varepsilon}$ as a volume form on the null surface one needs the covector $n_{\mu}$ to satisfy

[TABLE]

with $\varepsilon$ the volume 4-form, the one used to perform the integration in $\mathcal{V}$ in equation (3.57). This equation generalises the concept of orthonormal vector in the Gauss theorem and one can show that it implies eq. (3.57), see [32]. In our case, if one chooses to perform the integral over the null boundary as $\text{d}t\,\text{d}x\,\text{d}y$, which corresponds to a 3-form $\tilde{\varepsilon}_{\alpha\beta\gamma}$ which is completely antisymmetric in the variables $t,x$ and $y$, one has to normalise the orthogonal vector $n^{\mu}$ such that

[TABLE]

This is satisfied by the vector $n^{\mu}=(1,0,0,1)$ in the Minkowski coordinates $(t,x,y,z)$.

We now apply eq. (3.57) to the region depicted in Fig. 3 and use eq. (3.51). There is no dependence on $x$ and $y$, so we can factor out the surface $\text{d}x\,\text{d}y$:

[TABLE]

We dropped the contribution of the surface $\partial\mathcal{V}_{0}^{\prime}$ since all fields vanish on this surface. The LHS is the difference in energy between $t=T$ and $t=0$, while the RHS gives the flux of energy across $\partial\mathcal{V}_{0}$. As shown in Appendix C, neither $\bar{\gamma}$ nor $\Delta\gamma$ contribute to the energy flux across $\partial\mathcal{V}_{0}$ (intuitively GWs move parallel to this surface). Conversely, since $\pi$ only depends on $u$, it does not contribute to the difference in energy. Since the flux of energy is only due to $\pi$, which is constant on $\partial\mathcal{V}_{0}$, the flux is proportional to $T$. We see that the dependence $\Delta\gamma\propto T$ in the LHS is necessary for the cancellation. Notice that the sign of $\Delta\gamma_{ij}$ in (3.44) is the correct one: it implies that the amplitude of GW is decreasing. Indeed, since the total energy is just the sum of the kinetic energy of $\gamma$ and $\pi$, $\gamma$ must decrease in amplitude as $\pi$ grows. In Appendix C we check these statements and verify eq. (3.60).

4 Nonlinearities

We now want to look at the effects of $\pi$ non-linearities (non-linearities of $\gamma$ are suppressed by further powers of $M_{\rm Pl}$). In particular, we are going to study the effect of non-linearities on the exponential amplification of dark energy fluctuations that occurs in the narrow-resonance regime.

4.1 $m_{3}^{3}$-operator

We start from the operator (1.3). In this case, the cubic self-interaction in the Lagrangian is

[TABLE]

where $\Box\equiv-\partial_{t}^{2}+\nabla^{2}$ and

[TABLE]

At early times, when $\pi$ is in the vacuum state, it is safe to neglect these terms. However, as $\pi$ grows their importance increases and they become comparable to the resonance term

[TABLE]

Since $M\ll\Lambda$, we expect this to happen rather quickly. Comparing the above operators, this takes place when $\Box\pi\sim({M^{3}}/{\Lambda^{2}})\dot{\gamma}_{ij}\sim\sqrt{\alpha}{H}\dot{\gamma}_{ij}$, which can be written as

[TABLE]

by using $\dot{\gamma}_{ij}\sim\omega M_{\rm Pl}h_{0}^{+}$ and $\Box\pi\sim(\omega/c_{s})\partial_{i}\pi$. When this happens, both the energy density of $\pi$ and the modification of the GW, see eq. (3.34), are small. Indeed, using the above equation one finds

[TABLE]

After this point one can no longer trust the Mathieu equation and the solutions for $\pi$ used earlier, eq. (3.36): non-linear terms change the fundamental frequency of the oscillator in eq. (3.13), so that originally unstable modes are driven out of their instability bands and a more sophisticated analysis is required. The same conclusion can be obtained from the Boltzmann analysis of App. B. When the resonance term is modified by $\sim\mathcal{O}(1)$ corrections, particles are produced also outside the thin-shell of momenta $\Delta k$. This dispersion in momentum space implies that the Bose-enhancement factor, and in turn the growth index $\mu_{k}$, are affected.

The above estimate is also in agreement with numerical results in the preheating literature (see e.g. [33, 34]). These simulations suggest that even small self-interactions of the produced fields are enough to qualitatively change the development of the resonance. In these works it was also shown that attractive potentials for the reheated particles modify, and eventually shut down, the exponential growth found at linear level. This conclusion is not surprising, since in these cases large field expectation values contribute positively to their effective mass, making the decay kinetically disfavoured as soon as large field values are reached. This result cannot be applied to our case because the derivative self-interactions in eq. (4.1) do not enter with a definite sign in the action. To reach a definitive answer, in the narrow resonance regime, one would need a full numerical analysis.

4.2 $\tilde{m}_{4}^{2}$-operator

We will now take $m_{3}^{3}=0$ and focus on the self-interaction of $\pi$ generated by the operator of eq. (1.4). This choice results technically natural since $m_{3}^{3}$ corresponds, in the covariant theory, to cubic Horndeski operators that feature a weakly-broken Galilean invariance (for details, see the discussion in [35, 36, 1]).

Clearly, on sub-Hubble scales the most important nonlinearities are due to operators containing two derivatives per field. Since we are interested in the regime $\alpha_{\rm H}\ll 1$, in this case the most relevant non-linearities are not cubic but quartic. Indeed cubic non-linearities are suppressed by $\alpha_{\rm H}/\alpha^{3/2}\cdot\partial^{2}\pi/\Lambda_{3}^{3}$, while the quartic ones by $\alpha_{\rm H}/\alpha^{2}\cdot(\partial^{2}\pi/\Lambda_{3}^{3})^{2}$. Thus, for $\alpha_{\rm H}\ll 1$, quartic non-linearities become relevant for a smaller value of $\partial^{2}\pi/\Lambda_{3}^{3}$ (unless $\alpha$ is huge, we will not be interested in this regime below). Notice that the cut-off of the theory is thus of order $\Lambda_{3}\alpha^{1/3}\alpha_{\rm H}^{-1/6}$. This scale is much larger than $\omega$ for the values of $\alpha_{\rm H}$ we are interested in, so that the GW experimental results are well within the regime of validity of the theory. See Fig. 4 for a comparison of the scales in the problem.

Let us start with a naive estimate for $c_{s}\sim 1$. The resonance will be affected by non-linearities when

[TABLE]

Going through the same calculations as in the case of the $m_{3}^{3}$-operator one gets the constraint

[TABLE]

For typical values of the parameters the RHS of this inequality is at most of order unity. This would mean that non-linearities are important when the GW signal is substantially modified.

Actually, it turns out that the estimate above is not quite correct, as a consequence of the detailed structure of the quartic Galileon interaction. We want to evaluate the importance of the interactions on the modes that grow fastest, i.e. on the saddle point. Using $(u,s,x,y)$ coordinates, one has $p_{s}=0$ on the saddle, see eqs. (3.6) and (3.25). This can be understood from the equation of motion, eq. (3.11): for a given frequency of a $\pi$ wave, i.e. for a given momentum $|\boldsymbol{p}|$, one maximises the forcing term if $p_{s}=0$. Therefore in these coordinates the derivative with respect to $s$ vanishes. However, since the coordinate $u$ is null, the inverse metric satisfies $g^{uu}=0$. This means that also the derivative with respect to $u$ will not appear in the interaction (there are only cross terms $\propto\partial_{u}\partial_{s}$). Therefore we are left only with the two coordinates $x$ and $y$. However, the structure of the quartic Galileon is such that it vanishes if one has less than three dimensions. We conclude that the quartic $\pi$ self-interaction vanishes on the saddle point. To get a correct estimate one is forced to look at deviations from the saddle, i.e. one has to estimate what is the typical range of $p_{s}$ around $p_{s}=0$ that contributes to the integrals like eq. (3.40).

The rms value of the various variables can be read from the inverse of the Hessian matrix eq. (3.26). The momentum $p_{s}$ can be written in terms of the integration variables $\{\xi,\Omega,\varphi\}$ as

[TABLE]

One can write

[TABLE]

Evaluated at the saddle, the matrix of the first derivatives of $p_{s}$ is zero, except for the entry $(\Omega,\Omega)$ which is given by

[TABLE]

The entry $(\Omega,\Omega)$ of the inverse of the Hessian reads $\left(\text{H}^{-1}\right)^{\Omega\Omega}=\frac{(1-c_{s}^{2})^{2}}{4\beta\tau}(4-\beta^{2})^{3/2}$ . Therefore we obtain

[TABLE]

Using this estimate one can revise the bound above as

[TABLE]

It is important to stress that the coefficient of the quartic self-interaction of $\pi$ is tied by symmetry with the one of the operator $\ddot{\gamma}_{ij}\partial_{i}\pi\partial_{j}\pi$, so that the effect of self-interactions cannot be suppressed. This statement holds even considering models with $c_{\rm T}\neq 1$: since we are considering a regime of very small $\alpha_{\rm H}$, comparable values for $c_{\rm T}-1$ are not ruled out experimentally. If one does not impose the constraints on the speed of GWs, instead of the single operator of eq. (1.4) one has two independent coefficients

[TABLE]

In terms of these parameters, the coefficient of the quartic Galileon was calculated in [37] and it reads

[TABLE]

On the other hand the coefficient of the operator $\ddot{\gamma}_{ij}\partial_{i}\pi\partial_{j}\pi$ reads (see App. B of [1])

[TABLE]

The two coincide, modulo a factor of 2, up to relative corrections suppressed by $c_{\rm T}^{2}-1$.

5 Observational signatures for $\tilde{m}_{4}^{2}$

5.1 Fundamental frequency

In the narrow resonance regime, $\beta\ll 1$, with the replacement $\Lambda^{2}\rightarrow\Lambda_{\star}^{3}\omega^{-1}$, from eq. (3.44) one gets a relative modification of the GW given by

[TABLE]

The amplitude of GWs in the plus polarization is given by [31]

[TABLE]

where $r$ is the distance from the source and $M_{c}$ is the chirp mass of the binary system. The number of cycles of the gravitational wave is given by

[TABLE]

where we have defined the gravitational wave frequency $f\equiv\omega/2\pi$. In our calculations we took a time-independent frequency, while in reality $f$ increases with time, during the binary inspiralling. The frequency can be taken as roughly constant for the number of cycles $N_{\rm cyc}$ required for $f$ to double in size. In particular, the number of cycles between $f$ and $2f$ is given by [31]

[TABLE]

Inverting these relations we can express the chirp mass as a function of the number of cycles and the frequency. Using this in eq. (5.2) we obtain

[TABLE]

In the calculation we approximate the GW amplitude as a constant, therefore we have to limit $v\lesssim r$.

Expressing $\beta$ in term of $\alpha_{\rm H}$, replacing $h_{0}^{+}$ from the above relation and using $v=r$, eq. (5.1) becomes

[TABLE]

Note that this expression is independent of $\bar{N}_{\rm cyc}$. Sizeable effects in the GW waveform can be obtained when the argument of the exponential is $\sim{\cal O}(10^{2})$, which translates into

[TABLE]

We have also to impose the constraint of eq. (4.12), which using eq. (5.5) reads

[TABLE]

where the term $\beta\bar{N}_{\rm cyc}$ roughly coincides with the argument of the exponential in eqs. (5.1) and (5.6). Finally, a further constraint to impose is the narrow resonance condition $\beta\ll 1$.

In Figs. 5(a) and 5(b) we plot these three constraints as a function of $\alpha_{\rm H}$ and the distance between the source and the resonant decay of the GW into $\pi$ (therefore this is not the distance between the source and the detector). The first plot is done for a ground-based interferometer with LIGO/Virgo-like sensitivity; in order to maximise the number of oscillation we choose a neutron-star event similar to GW170817 [7]. The second plot is for a space-based interferometer with LISA-like sensitivity [38] and a binary black hole event similar to GW150914 [39]. The blue region corresponds to a sizeable modification of the GW signal, calculable within our approximation. The neutron-star merger GW170817 is at a distance of approximately 40 Mpc. The absence of sizeable effects ($\Delta\gamma>0.1\gamma$; of course future measurements will improve this sensitivity) in the observed event puts constraints on a resonant effect that takes place at less than 40 Mpc from the source and rules out the interval: $3\times 10^{-20}\lesssim\alpha_{\rm H}\lesssim 10^{-18}$. Future measurements by LISA of events similar to GW150914 can, on the other hand, constrain the range $10^{-16}\lesssim\alpha_{\rm H}\lesssim 10^{-10}$.

It may be useful to summarize here the set of assumptions for the validity of the constraints in these figures:

• •

Narrow resonance regime, i.e. $\beta\sim\alpha_{\rm H}(\omega/H)^{2}h_{0}^{+}\lesssim 1$. This is diplayed by the blue-dashed line in the figures.

• •

The modification of the GW is sufficiently sizable to be observed, i.e. $|\Delta\gamma|>0.1\gamma$. This is diplayed by the purple-dashed line in the figures. The modification takes place outside the Vainshtein radius, which for the values of $\alpha_{\rm H}$ that we are considering corresponds to distances from the source smaller than those shown on the horizontal axis. Moreover, to be relevant the modification must occur before the detection.

• •

Nonlinearities in $\pi$ are small so that their effect on the Mathieu equation is negligible, see discussion in Sec. 4.2. This is displayed by the red-dashed lines in the figures.

• •

In the figures we assume $c_{s}=1/2$ and the constraints would change for different values of $c_{s}$. From eq. (3.44) the effect is proportional to $(1-c_{s}^{2})^{2}$ so that it is suppressed for values of $c_{s}$ close to 1. Moreover, as we discussed in Sec. 3.1, for $c_{s}-1\lesssim 10^{-14}$ the $\pi$ lightcone is too wide and our calculations do not apply. (Notice also that in our estimates of the various scales we are assuming that $c_{s}$ is not parametrically small.)

For comparison, the bound coming from the perturbative decay of the graviton, see eq. (46) of [1], reads $\alpha_{\rm H}\lesssim 10^{-10}$. Such values of $\alpha_{\rm H}$ are outside the narrow resonance regime studied in this paper. One expects that the effect of large occupation number and nonlinearities of $\pi$ are also important for larger values of $\alpha_{\rm H}$ so that this perturbative bound should be revised. Moreover, notice that when the GW is closer to the source, its amplitude is larger and $\beta$ will exceed unity. At a certain point one enters the Vainshtein regime and the coupling $\gamma\pi\pi$ is reduced (since we are considering very small values of the coupling, the Vainshtein radius will be correspondingly small). We will study all these aspects in a forthcoming publication [13].

5.2 Higher harmonics and precursors

Since for the operator $\tilde{m}_{4}^{2}$ one can trust the parametric growth and the change in the GW signal $\Delta\gamma$, we want to study this correction in more detail. At leading order in $q$ the change $\Delta\gamma$ has the same frequency as the original wave. However, corrections to the leading solution (3.36) introduce higher harmonics in the GW signal. Higher harmonics appear in two quantitatively different ways. First, we can consider higher instability bands, see Fig. 2(a). However, in the narrow resonance approximation the instability coefficient in the $m$-th band scales as $\mu_{m}\sim q^{m}/(m!)^{2}\sim\beta^{m}$ [29, 30]. Therefore, in the long $\tau$ limit they are exponentially suppressed compared to the fundamental band.

Second, we can stay in the first band and look at the modification of the $\pi$ solution at higher orders in $q$. This will give an effect which is only suppressed by positive powers of $q$. If we include the first two corrections to eq. (3.36), we have666In this section we write the formulas for the case of $m_{3}^{3}$, so that they can be easily compared with Section 3.3. Notice that, since in this case $\pi$ is coupled with $\dot{\gamma}$, while in the case $\tilde{m}_{4}^{2}$ it is coupled with $\ddot{\gamma}$, there is an overall shift of $\pi/2$ between the two cases.

[TABLE]

Clearly, in this case the higher frequency component grows as fast as the original frequency. Therefore it is not necessary to specify the correction $\delta\mu$. The correction to $|f_{\tilde{p}}|^{2}$ at order $q^{2}$ is thus

[TABLE]

On the other hand, when $q$ becomes negative the solution can be obtained by performing the transformation $\tau\rightarrow\tau+\pi/2$, or equivalently $q\rightarrow-q$ together with $\sigma\rightarrow-\sigma-\pi/2$. By using the latter transformation and by recalling that by flipping the sign of $q$ we send the decreasing mode into the growing one, we get the contribution for negative $q$:

[TABLE]

By direct calculation we can also show that, as in (3.41), $|c_{+}|=|c_{+}^{\prime}|=1$ at lowest order in $q$. We also note that, since all the frequencies in the expressions above are multiples of $\omega$, the position of the saddle point is not affected and it is simply fixed by the exponential function. Hence the expression for $J(u)$ can be obtained by using the saddle-point procedure as in (3.43). We obtain

[TABLE]

(The result for the case of $\tilde{m}_{4}^{2}$, which is the focus of this section, can be obtained by the replacement $\Lambda^{2}\to\Lambda_{\star}^{3}\omega^{-1}$.) We are interested in the $\tau$-dependent oscillatory part of (5.12). Compared to (3.43) we get a correction at the fundamental frequency linear in $\beta$ (with no phase-shift). Moreover, to find the first higher harmonics we have to go to order $\beta^{2}$, where indeed we encounter a term with three-times the frequency $\omega$. Note that by going to the second band instead, we would find higher harmonics of frequency $2\omega$.

Higher harmonics enter in the bandwidth of the detector at earlier time compared to the main signal. In some sense they are precursors of the main signal. To assess their potential observability, one has to consider that also the post-Newtonian expansion generates terms with a frequency which is a multiple of the original one. We leave this study to future work.

6 Conclusions and future directions

We have studied the effects of a classical GW on scalar field fluctuations, in the context of dark energy models parametrized by the EFT of Dark Energy. The GW acts as a classical background and modifies the dynamics of dark energy perturbations, described by $\pi$, leading to parametric resonant production of $\pi$ fluctuations. This regime is described by a Mathieau equation, where modes in the unstable bands grow exponentially. In the regime of narrow resonance, corresponding to a small amplitude of the GW, the instability can be studied analytically in detail. One can also include the back-reaction of the produced $\pi$ on the original GW, which leads to a change in the signal. The resonant growth is however very sensitive to the non-linearities of the produced fields: when non-linearities of $\pi$ are sizeable the resonance is damped. This regime cannot be captured analytically: numerical simulations are needed to confirm this effect. This happens for the operator $m_{3}^{3}$, while a sizeable modification of the GW signal is possible for the operator $\tilde{m}_{4}^{2}$, at least in some range of parameters, see Figs. 5(a) and 5(b).

Many questions remain open and will be the subject of further studies. For larger amplitudes of the GW, one exits the regime of narrow resonance: the GW may induce instabilities, which will be studied in [13]. This will also allow to revise the robustness of the perturbative calculation of [1]. In the same paper we will also study how these effects are changed in the presence of Vainshtein screening. Simulations are probably required if one wants to study the regime of parametric resonance in the presence of sizeable $\pi$ self-interactions. Similar effects are at play in preheating after inflation and are also addressed numerically, see e.g. [33, 34]. It would also be interesting to try to envisage methods to look for the small changes in the GW signal. For instance, the generation of higher harmonics that enter the detection bandwidth before the main signal is a striking possible effect if it can be disentangled from post-Newtonian corrections. The effects that we studied look quite generic to all theories in which gravity is modified with a quite low cut-off. Therefore, it would be nice to explore other setups, starting from the DGP model. This model features the same non-linear interaction of $\pi$ as the ones of the $m_{3}^{3}$ operator; however the extra-dimensional origin changes the dynamics of GWs so our analysis cannot be applied directly. We hope to come back to all these points in the future.

Acknowledgements

It is a pleasure to thank M. Celoria, G. Cusin, P. Ferreira, M. Lagos, M. Lewandowski, S. Melville, A. Nicolis, and F. Schmidt for useful discussions. We also thank the anonymous referee for a very careful reading of the article and insightful comments. P. C. and F. V. acknowledge the kind hospitality of PUCV, Valparaiso, Chile where part of this work was carried on.

Appendix A Interactions in Spatially-Flat Gauge

In this appendix we are interested in redoing the computation of Section 2.1 in the spatial-flat gauge. Here the metric is written in the general decomposition (see e.g. [40])

[TABLE]

with $\tilde{\gamma}_{ij}$ being transverse, $\tilde{\partial}_{i}\tilde{\gamma}_{ij}=0$, and traceless $\delta_{ij}\tilde{\gamma}_{ij}=0$. The shift vectors $\tilde{N}_{i}$ can be decomposed as $\tilde{N}_{i}=\tilde{\partial}_{i}\tilde{\psi}+\hat{\tilde{N}}_{i}$ where $\tilde{\partial}_{i}\hat{\tilde{N}}_{i}=0$ ($\tilde{\partial}_{\mu}\equiv\partial/\partial\tilde{x}^{\mu}$). Note that in this gauge $\tilde{\pi}(\tilde{x})$ denotes the Goldstone field. Since both $\delta N$ and $\tilde{N}_{i}$ are Lagrange multipliers in this spatially-flat gauge, by solving the constraint equations (see e.g. [1, 41]) they can be expressed up to first order as

[TABLE]

and

[TABLE]

where $\tilde{\nabla}^{2}\equiv\tilde{\partial}_{i}\tilde{\partial}_{i}$. Notice that we have kept the non-local term in (A.3) since this will contribute to the vertex that we are going to consider. According to [1], $\tilde{\pi}$ is canonically normalised as

[TABLE]

while the canonical normalisation of $\tilde{\gamma}_{ij}$ remains the same as in the Newtonian gauge. Here we keep a distinction between the canonical field $\pi_{c}$ and the non-canonical field $\tilde{\pi}$ unlike in the main text.

We are going to investigate the non-linear $\gamma\pi\pi$ vertex which contains three or more derivatives (we will drop all the rest). In principle, this vertex gets contributions from the action (1.2) through the Einstein-Hilbert and the $S_{m_{3}}$ terms. Let us first consider the contribution from $S_{m_{3}}$ term. As usual, the extrinsic curvature is defined by

[TABLE]

where $\tilde{D}_{i}$ denotes the 3d covariant derivative with respect to the induced metric $\tilde{h}_{ij}$. It is straightforward to show that a variation of $\tilde{K}$ ($\tilde{K}=\tilde{h}^{ij}\tilde{K}_{ij}$) subjected to the metric (A.1) contains $a^{-2}\tilde{\gamma}_{ij}\tilde{\partial}_{i}\tilde{\partial}_{i}\tilde{\psi}$. After performing the Stuekelberg trick of eq. (2.1) and (2.2), the contribution from $S_{m_{3}}$ which has three derivatives reads

[TABLE]

where we have taken the term $-2(1-\alpha_{N})\dot{\tilde{\pi}}$ from $\delta\tilde{g}^{00}$ and used (A.2)-(A.3).

Using the canonical normalisations both for $\tilde{\pi}$ (A.4) and $\tilde{\gamma}_{ij}$ one obtains the same interaction that we have obtained in the Newtonian gauge (2.8). Thus we expect that the contribution arising from $S_{EH}$ cancels out. To show this, notice that the contribution from the Einstein-Hilbert term comes from $N(\tilde{K}_{ij}\tilde{K}^{ij}-\tilde{K}^{2})$. More specifically, we have

[TABLE]

where we have performed a few integration by parts. Using (A.2) and (A.3) in $S_{\text{EH}}$ one obtains

[TABLE]

where we have used the linear equations of motion for $\tilde{\gamma}_{ij}$, $\ddot{\tilde{\gamma}}_{ij}+3H\dot{\tilde{\gamma}}_{ij}-\frac{1}{a^{2}}\tilde{\nabla}^{2}\tilde{\gamma}_{ij}=0$, and for $\tilde{\pi}$, $\ddot{\tilde{\pi}}+3H\dot{\tilde{\pi}}-\frac{c_{s}^{2}}{a^{2}}\tilde{\nabla}^{2}\tilde{\pi}=0$.

Notice that the first two terms on RHS vanish due to the fact that $\alpha_{N}+\alpha_{\psi}=0$. Let us consider the prefactor of the last term. Using the expression of $c_{s}^{2}$ (2.13) and the definitions of all those $\alpha$-parameters (A.2)-(A.3), the prefactor can be rewritten as

[TABLE]

In our case the coupling with matter has been neglected ($\rho_{m}=P_{m}=0$), therefore the parameter c is equal to $-M^{2}_{Pl}\dot{H}$ (see e.g. [3, 5]). The prefactor is then given by

[TABLE]

which is again zero since $\alpha_{N}+\alpha_{\psi}=0$ in this case. Therefore $S_{EH}$ gives no contribution to our $\gamma\pi\pi$ vertex and $S_{m_{3}}$, on the other hand, gives the same result that we have obtained in the Newtonian gauge (2.8).

Appendix B Parametric resonance as Bose enhancement

In this appendix we want to reinterpret the exponential growth due to parametric resonance as the Bose enhancement of the perturbative decay $\gamma\rightarrow\pi\pi$. To see this, we study the Boltzmann equation for the number density of dark energy fluctuations. We denote by $n_{\boldsymbol{k}}^{\pi}$ and $n_{\boldsymbol{k}}^{\gamma}$ the occupation numbers, respectively, of $\pi$ and $\gamma$. Moreover, the number density for the particle species $\pi$ is defined as

[TABLE]

and an analogous definition holds for $\gamma$.

Let us consider a collection of gravitons with frequency $\omega$, each of them decaying into two $\pi$-particles. For concreteness, we will focus on the case $c_{s}^{2}\ll 1$, for which the two momenta of $\pi$, $\boldsymbol{k}$ and $-\boldsymbol{k}$, have opposite directions and equal magnitudes $k=\omega/(2c_{s})$. Following [42] and denoting by $\Gamma_{\gamma\rightarrow\pi\pi}$ the tree-level decay rate (see e.g. (2.14)), the rate of change of $n_{\pi}$ in a given volume $V$ is

[TABLE]

where the factor of 2 accounts for two identical particles in the final state. As explained earlier in Sec. 3.1, for small $c_{s}$ we can use $u$ as time. On the right-hand side, we have neglected integration over the angle $\varphi$ which would appear when considering only one of the GW polarizations. For $n_{\boldsymbol{k}}^{\pi}=n_{-\boldsymbol{k}}^{\pi}=n_{k}^{\pi}$ and $n_{\omega}^{\gamma}\gg\{n_{\boldsymbol{k}}^{\pi},1\}$ we find

[TABLE]

where we have introduced the number density of gravitons, here given by $n_{\gamma}=n_{\omega}^{\gamma}/V$.

The produced $\pi$-particles end up populating the spherical shell $k=k_{0}\pm\Delta k/2$ of radius $k_{0}\simeq\omega/(2c_{s})$ and thickness $\Delta k$, so that their occupation number is related to their number density by

[TABLE]

The thickness is given by comparing the time-independent part of the equation of motion for $\pi$, eq. (3.3), with the amplitude of its periodic part. Using $\Omega=c_{s}\ll 1$ valid for small $c_{s}$, we obtain

[TABLE]

Plugging the expressions of $k_{0}$ and $\Delta k$ into (B.4) and using $n_{\gamma}\simeq\omega(M_{\rm Pl}h^{+}_{0})^{2}$, the occupation number can be written as

[TABLE]

where we have focused on the case of the operator $m_{3}^{3}$ (which can be straightforwardly extended to the case of the operator $\tilde{m}_{4}^{2}$ by the replacement $\Lambda^{2}\to\Lambda_{\star}^{3}\omega^{-1}$, see discussion in Sec. 2.2) and used the definition of $\beta$, eq. (2.18), to rewrite $(M_{\rm Pl}h^{+}_{0})^{2}$. The Bose condensation effect becomes important for $n_{k}^{\pi}\gg 1$ or, using the above equation, for

[TABLE]

In this case, we can solve eq. (B.3) with the decay rate given by eq. (2.14). This gives

[TABLE]

which displays an instability similar to that encountered above in eq. (3.27), but with a different exponent. Notice that the approach of Appendix B is approximate and does not reproduce the correct numerical factors in the timescale of the instability. Of course, a more precise calculation would give the same answer.

It is useful to check that our formula for the modification of the GW, eq. (3.44), smoothly interpolates with the perturbative decay result, eq. (2.14), when the occupation number becomes small. In the regime in which Bose enhancement is negligible $n_{\pi}/n_{\gamma}\sim\Gamma u$, see Fig. 1. Using this in eq. (B.6) we get $n^{\pi}_{k}\sim\beta\omega u$. Not surprisingly this is the parameter that enters the exponential growth of the instability. Our saddle-point treatment is valid for $\beta\omega u\gg 1$, but we expect that, when $\beta\omega u\sim 1$, it gives a result of the same order as the perturbative decay of gravitons. Indeed if we plug this equality in eq. (3.44) one gets (for $c_{s}\ll 1$)

[TABLE]

This is indeed the perturbative result: the original GW, $M_{\rm Pl}h_{0}^{+}$, changes with a rate $\Gamma$ for a time of order $v$.

Appendix C Details on the conservation of energy

C.1 $m_{3}^{3}$-operator

In this Appendix we check the conservation of energy discussed in Section 3.5. First of all, let us verify that GWs do not contribute to the flux of energy across $\partial\mathcal{V}_{0}$ (see Fig. 3). This is clearly true for $\bar{\gamma}_{ij}$ but it holds at order $\Delta\gamma$ too. Indeed we have, using $\dot{\bar{\gamma}}_{ij}=-\partial_{z}\bar{\gamma}_{ij}$ and $\partial_{k}\bar{\gamma}_{ij}=\delta_{kz}\partial_{z}\bar{\gamma}_{ij}$,

[TABLE]

Let us now calculate the LHS of eq. (3.60), i.e. the variation of the total energy in the region. This is only due to $\Delta\gamma_{ij}$, since $\pi$ and $\bar{\gamma}_{ij}$ depend on $u$ only. One has

[TABLE]

The integral over $\partial\mathcal{V}_{2}\cup\partial\mathcal{V}_{1}$ reduces to

[TABLE]

The RHS of eq. (3.60) gets contribution only from $\pi$:

[TABLE]

The equation for energy conservation, eq. (3.60), becomes

[TABLE]

Note that the linear $v$ dependence of $\Delta\gamma$ is essential: if it was not the case then the two terms in (C.6) would have a different $T$ dependence, with no chance of adding up to zero. To explicitly check energy conservation one should integrate (C.6). A faster way to check this is to take a derivative with respect to $|z_{0}|$ (or equivalently $u$). The resulting equation can be shown to be satisfied by the solution for $\pi$, (3.13). After simplifying $T$ and taking the derivative with respect to $z_{0}$, (C.6) becomes

[TABLE]

The expression of $\langle T^{00}_{\pi}\rangle$ is given by eq. (3.2) while

[TABLE]

Therefore we can write

[TABLE]

The remaining term contains

[TABLE]

At this point by adding up these contributions we can collect the terms with ${\partial_{u}f_{\tilde{p}}^{\star}}$. They are

[TABLE]

This equation is solved by $f_{\tilde{p}}$: this is easily seen by comparing with equation (3.13).

C.2 $\tilde{m}_{4}^{2}$-operator

In this section we are going to use the same logic as the one in the previous section to verify the conservation of energy for $\tilde{m}_{4}^{2}$-operator. According to the Lagrangian (2.19), the components of the energy-momentum tensor are given by

[TABLE]

Notice that there is an extra term in $T^{0}_{\ 0}$ due to the interaction $\gamma\pi\pi$, unlike the case of $m_{3}^{3}$-operator. Since this new piece is second order in $\pi$, it can be approximated as $-\frac{2}{\Lambda_{\star}^{3}}\dot{\bar{\gamma}}_{ij}\partial_{i}\dot{\pi}\partial_{j}\pi$.

Let us first consider the RHS of (3.60), taking into account that $\dot{\bar{\gamma}}_{ij}=-\partial_{z}\bar{\gamma}_{ij}$ and $\partial_{k}\bar{\gamma}_{ij}=\delta_{kz}\partial_{z}\bar{\gamma}_{ij}$,

[TABLE]

where all terms involving only GWs added up to zero because of eq. (C.1). The LHS of (3.60) gets contributions only from $\Delta\gamma_{ij}$, as for the ${m}_{3}^{3}$-operator case. We have

[TABLE]

Integrating $T^{00}$ over $\partial\mathcal{V}_{2}\cup\partial\mathcal{V}_{1}$ gives

[TABLE]

similarly to eq. (C.3). Using (C.2) and (C.18), eq. (3.60) becomes

[TABLE]

Like in the previous section, one can verify this equation by taking a derivative with respect to $|z_{0}|$ (or equivalently $u$):

[TABLE]

As we have shown before, the first term of LHS can be expressed as

[TABLE]

Similarly, the second term can be rewritten as

[TABLE]

Adding up (C.21) and (C.22) together one therefore obtains

[TABLE]

This coincides with the equation of motion for $f_{\tilde{p}}$, which can be obtained by expanding eq. (2.23) in Fourier modes.