# Gravitational wave turbulence: a multiple time scale approach for quartic
wave interactions

Benoît Gay<EMAIL_ADDRESS>Sébastien Galtier
<EMAIL_ADDRESS>Laboratoire de Physique des Plasmas,
Université Paris-Saclay, CNRS, École Polytechnique, Sorbonne Université,
Observatoire de Paris, F-91128 Palaiseau, France

###### Abstract

Wave turbulence is by nature a multiple time scale problem for which there is
a natural asymptotic closure. The main result of this analytical theory is the
kinetic equation that describes the long-time statistical behaviour of such
turbulence composed of a set of weakly nonlinear interacting waves. In the
case of gravitational waves, it involves four-wave interactions and two
invariants, energy and wave action. Although the kinetic equation of
gravitational wave turbulence has been published with the Hadad-Zakharov
metric, along with their physical properties, the detailed derivation has not
been shown. Following the seminal work of Newell [1] for gravity/surface
waves, we present the multiple time scale method, rarely used to derive the
kinetic equations, and clarify the underlying assumptions and methodology.
This formalism is applied to a wave amplitude equation obtained using an
Eulerian approach. It leads to a kinetic equation slightly different from the
one originally published, with a wave equation obtained using a Hamiltonian
approach; we verify, however, that the two formulations are fully compatible
when the number of symmetries used is the same. We also show that the exact
solutions (Kolmogorov-Zakharov spectra) exhibit the same power laws and
cascade directions. Furthermore, the use of the multiple time scale method
reveals that the system retains the memory of the initially condition up to a
certain level (second order) of development in time.

Gravitational waves, Weak turbulence, Wave-turbulence interactions,
Homogeneous turbulence

##  I Introduction

The initial detection of gravitational waves by the LIGO collaboration [2]
came nearly a century after their original theoretical prediction by Albert
Einstein [3]. The observed deformations on Earth are of very small magnitude
($10^{-21}$), allowing for a predominantly linear approach to study their
propagation quickly after they were formed [4]. However, when examining the
history of the early universe, certain events may have given rise to
gravitational waves of higher amplitude. For instance, we can cite the first-
order phase transition [5, 6, 7] or the self-ordering of scalar fields [8].
Moreover, gravitational waves of significant amplitude are anticipated to have
been generated during cosmic inflation [9, 10], and numerous endeavours are
underway to detect them [11, 12]. The recent publication of the NANOGrav’s
survey [13] provides solid evidence for the existence of a gravitational-wave
background. While this result is consistent with a population of supermassive
black-holes, we cannot exclude the primordial cosmological origin of this
phenomenon.

When dealing with gravitational waves of large amplitude, the linear approach
becomes inappropriate and the introduction of non-linear effects inevitable.
Among the possible methods, wave turbulence offers a valuable set of tools for
examining the weakly non-linear regime, where the injection and dissipation of
an invariant (typically the energy or the wave action) take place across
different scales. The transfer of this invariant from a larger scale to a
smaller scale is referred to as a direct cascade, while the transfer from a
smaller scale to a larger scale is known as an inverse cascade. Such cascades
manifest in a diverse range of natural phenomena, including surface waves
[14], rotating fluids [15], plasma physics [16], vibrating plates [17], Bose-
Einstein condensate [18], optical waves [19] and gravitational waves [20].

The foundations of the turbulence theory were laid in particular by
Kraichnan’s direct interaction approximation, which elucidates the
conservation of energy in hydrodynamics through triadic wave interactions
[21]. The main developments in wave turbulence came a little later with Benney
and Saffman [22] who introduced a multiple time scale method [23] for a
problem involving three-wave interactions, and Benney and Newell [24, 1] who
extended the theory to quartic interactions with an application to gravity
(surface) waves. In both cases, the Eulerian approach was used and the main
result was the derivation of the kinetic equation that describes the temporal
evolution of spectral invariants such as energy. Simultaneously, the Soviet
school focused on wave turbulence, particularly in plasmas, using a
Hamiltonian approach. Employing a distinct method known as the random phase
approximation, they derived the kinetic equations [25, 26, 27, 28, 29], but
also found the finite flux spectra (called Kolmogorov-Zakharov spectra) as
exact stationary solutions [14, 30].

Over the past ten years, progress has been made in the field of non-linear
perturbative approaches to general relativity with a significant focus on the
anti-de Sitter space. In particular, it has been shown that this space is
unstable under arbitrarily small generic perturbations, leading to a diffusion
of energy towards higher frequencies and potentially to the formation of a
black hole [31]. Numerical investigations have confirmed the plausible
turbulent nature of spherically symmetric gravitational collapse in
asymptotically anti-de Sitter space, with a scalar curvature showing a power
law spectrum close to $-5/3$ [32]. This result was improperly called
Kolmogorov or Kolmogorov-Zakharov spectrum, whereas in turbulence each problem
has its own power law energy spectrum index ($-5/3$ being the scaling mainly
found in hydrodynamics), and the denomination Kolmogorov-Zakharov is only
given to a finite flux spectrum (exact) solution of the kinetic equation [33,
20, 34]. Note that the conclusion on the stability study is not universally
applicable as other solutions can remain non-linearly stable [35, 36].
Interestingly, it is for such an anti-de Sitter space that a technique
inspired by the renormalization group theory has been introduced to deal with
secular terms that naturally appear in a perturbation theory [37]. This
technique for treating secular terms bears some resemblance to the multiple
time scale method that we will present in this paper. Other investigations
[38] have revealed the existence of conserved quantities such as energy or
wave action, the latter not well known to non-specialists but often found in
wave turbulence. For example, in gravitational wave turbulence, these are two
conserved quantities that offer further insights into the dynamics [20].

It is only recently that significant progress has been made in understanding
gravitational wave turbulence, primarily spearheaded by Galtier and Nazarenko.
After proving that this turbulence necessarily occurs through quartets of
interactions, they used a Hamiltonian approach based on the metric proposed by
Hadad and Zakharov [39], which provides a concise and consistent set of
equations, to derive the kinetic equation and its exact statistical solutions.
They identified two distinct cascades: a direct cascade in energy and an
inverse cascade in wave action [20]. The latter turned out to be an explosive
cascade able in principle of reaching the wave number $\mathbf{k}=\mathbf{0}$
in a finite time [40]. However, the wave turbulence description fails at scale
$k_{s}$ where turbulence becomes strongly non-linear. Using a phenomenological
model, it is found that the wave action spectrum undergoes a change in slope
that remains compatible with an explosive inverse cascade. In the presence of
a small scale excitation, a condensate is expected to form rapidly at small
wave numbers. This scenario was then proposed to model the primordial
universe: the idea is that the growth of the condensate can be linked to the
expansion factor of the Robertson-Walker model [41]. The explosive growth
would then correspond to a phase of inflation. The associated dilution, and
the decrease of the nonlinearities, provides a natural mechanism for stopping
the expansion, leaving the universe with a fossil spectrum corresponding to
the Harrison-Zeldovich spectrum which is compatible (at main order) with the
observations [42]. As the contribution of energy tends to zero in the small
$k$ limit, this turbulent inflation scenario therefore requires negligible
energy. More recently, using direct numerical simulations, the first direct
evidence of a dual cascade has been found [43], proving the relevance of the
analytical theory of gravitational wave turbulence.

The aim of this paper is to establish the kinetic equation for gravitational
wave turbulence using a multiple time scale approach initially proposed for
three-wave interactions [22] and then generalized to quartic interactions [24,
1]. This demonstration was not given in the original paper [20]. As turbulence
is an obscure subject in cosmology, we think it is important to explain in a
pedagogical way (the derivation is long and non-trivial) how to derive the
kinetic equation. Section II will serve as a reminder of the preliminary
assumptions required to derive the evolution equation for the fields within
the Hadad-Zakharov metric [39]. Here, we follow an Eulerian approach, which is
new. Subsequently, in section III, we will describe the general process of
deriving the kinetic equation. This a technical section in which many details
are given. In section IV, we will focus on the specific case of gravitational
waves, whose symmetries will have been previously identified. Finally, a
conclusion is proposed in section V.

##  II Expansion of the Einstein equation

Throughout the entire paper, we adopt a unit system where the speed of light,
$c$, is set to 1. Consequently, the dispersion relation for gravitational
waves takes the simple form: $\omega_{\mathbf{k}}=k$, where $k$ represents the
magnitude of the wave vector $\mathbf{k}$. While we do not employ the Einstein
summation convention, we use the following convention:

$\delta^{\alpha_{1}\dots\alpha_{p}}_{\beta_{1}\dots\beta_{q}}(\mathbf{k})=\delta\left(\sum_{i=1}^{p}\mathbf{k}_{\alpha_{i}}-\sum_{i=1}^{q}\mathbf{k}_{\beta_{i}}\right),$
(1)

and

$\Omega^{\alpha_{1}\dots\alpha_{p}}_{\beta_{1}\dots\beta_{q}}=\sum_{i=1}^{p}s_{\alpha_{i}}\omega_{\mathbf{k}_{\alpha_{i}}}-\sum_{i=1}^{q}s_{\beta_{i}}\omega_{\mathbf{k}_{\beta_{i}}},$
(2)

where $s=\pm 1$ is the directional polarity and
$\left(\alpha_{i}\right)_{1\leq i\leq p}$ and $\left(\beta_{i}\right)_{1\leq
i\leq q}$ are two family of indices.

This notation allows for concise representation and facilitate the clarity of
the discussion. Below, we provide an overview of the wave amplitude
calculation as outlined in [20].

###  II.1 Hadad–Zakharov metric

We begin by considering an observer located in the vacuum at a significant
distance from the system under study. This observer utilizes a Cartesian
coordinate system $(x,y,z)\in\mathbb{R}^{3}$ and its own proper time
$t\in\mathbb{R}_{+}$ to describe the interactions. Following [39], we assume
that the metric is diagonal and independent of the $z$ coordinate. Thus,
$\partial_{z}$ is a Killing field and the metric can be re-written as the
following:

$g_{\mu\nu}=e^{-2\varphi}\begin{pmatrix}-(1+\gamma)^{2}&0&0&0\\\
0&(1+\beta)^{2}&0&0\\\ 0&0&(1+\alpha)^{2}&0\\\ 0&0&0&e^{4\varphi}\\\
\end{pmatrix},$ (3)

where $\alpha,\beta,\gamma$, and $\varphi$ are real functions that depend
solely on $t$, $x$, and $y$. It is important to note that this metric does not
take into account the two polarizations of gravitational waves. This
simplification considerably reduces the complexity of the system under study,
while retaining the main ingredients of wave turbulence: waves and
nonlinearities.

The vacuum Einstein equation $R_{\mu\nu}=0$, with $R_{\mu\nu}$ the Ricci
tensor of the metric $g_{\mu\nu}$, yields a set of ten equations. However, as
demonstrated in [39], three of these equations are trivially zero, and three
are redundant with the remaining four equations. Consequently, we can focus on
the following four non-linear equations:

$\displaystyle\partial_{x}\partial_{t}\alpha=-$ $\displaystyle
2(1+\alpha)\partial_{t}\varphi\partial_{x}\varphi+\frac{\partial_{t}\beta}{1+\beta}\partial_{x}\alpha+\frac{\partial_{x}\gamma}{1+\gamma}\partial_{t}\alpha,$
(4a) $\displaystyle\partial_{y}\partial_{t}\beta=-$ $\displaystyle
2(1+\beta)\partial_{t}\varphi\partial_{y}\varphi+\frac{\partial_{t}\alpha}{1+\alpha}\partial_{y}\beta+\frac{\partial_{y}\gamma}{1+\gamma}\partial_{t}\beta,$
(4b) $\displaystyle\partial_{x}\partial_{y}\gamma=-$ $\displaystyle
2(1+\gamma)\partial_{x}\varphi\partial_{y}\varphi+\frac{\partial_{x}\alpha}{1+\alpha}\partial_{y}\gamma+\frac{\partial_{y}\beta}{1+\beta}\partial_{x}\gamma,$
(4c) and
$\partial_{t}\left[\frac{(1+\alpha)(1+\beta)}{1+\gamma}\partial_{t}{\varphi}\right]=\partial_{x}\left[\frac{(1+\alpha)(1+\gamma)}{(1+\beta)}\partial_{x}\varphi\right]+\partial_{y}\left[\frac{(1+\beta)(1+\gamma)}{(1+\alpha)}\partial_{y}\varphi\right].$
(4d)

The first three expressions correspond to three constraint equations, while
the last one describes the non-linear dynamics of the metric.

In the linear case, the fields $\alpha$, $\beta$, and $\gamma$ are assumed to
be null (this is also our initial condition at $t=0$) and consequently only
the following two-dimensional wave equation remains:

$\partial_{t}\partial_{t}\varphi-\partial_{x}\partial_{x}\varphi-\partial_{y}\partial_{y}\varphi=0.$
(5)

The solutions of (5) are expressed as a superposition of plane waves:

$\varphi(\mathbf{x},t)=\int_{\mathbb{R}^{2}}\left[A^{+}(\mathbf{k})e^{-i(\omega_{\mathbf{k}}t-\mathbf{k}\cdot\mathbf{x})}+A^{-}(\mathbf{k})e^{-i(\omega_{\mathbf{k}}t+\mathbf{k}\cdot\mathbf{x})}\right]~{}\mathrm{d}^{2}\mathbf{k},$
(6)

where by definition $\mathbf{x}=\left(x,y\right)$,
$\mathbf{k}=\left(p,q\right)$ and $A^{\pm}(\mathbf{k})$ are square-integrable
functions mapping $\mathbb{R}^{2}$ to $\mathbb{R}^{2}$.

To focus on small oscillations, we introduce a small parameter $\epsilon$
(with $0<\epsilon\ll 1$) which measures the amplitude of the field $\varphi$.
The $\alpha$, $\beta$, and $\gamma$ fields can be expanded around their
equilibrium positions. At first nonlinear order, equations (4d) can be
simplified as follows:

$\displaystyle\partial_{x}\partial_{t}\alpha$
$\displaystyle=-2\partial_{t}\varphi\partial_{x}\varphi,$ (7a)
$\displaystyle\partial_{y}\partial_{t}\beta$
$\displaystyle=-2\partial_{t}\varphi\partial_{y}\varphi,$ (7b)
$\displaystyle\partial_{x}\partial_{y}\gamma$
$\displaystyle=-2\partial_{x}\varphi\partial_{y}\varphi,$ (7c) and
$\partial_{t}\left[\left(1+{\alpha}+{\beta}-{\gamma}\right)\partial_{t}\varphi\right]=\partial_{x}\left[\left(1+{\alpha}-{\beta}+{\gamma}\right)\partial_{x}\varphi\right]+\partial_{y}\left[\left(1-{\alpha}+{\beta}+{\gamma}\right)\partial_{y}\varphi\right].$
(7d)

At this stage, the constraint equations (7a–7c) clearly show that the
$\alpha$, $\beta$, and $\gamma$ fields depend quadratically on $\varphi$.
Therefore, the non-linearities involved in the dynamic equation (7d) are cubic
in nature so that we have to deal with four-wave interactions in Fourier
space.

Note that we consider a continuous medium which can lead to mathematical
difficulties connected with infinite dimensional phase spaces. For this
reason, it is preferable to assume a variable spatially periodic over a box of
finite size $L$. However, in the derivation of the kinetic equation, the limit
$L\to+\infty$ is finally taken (before the long time limit, or equivalently
the limit $\epsilon\to 0$). As both approach leads to the same kinetic
equation, for simplicity, we anticipate this result and follow the original
approach of [22].

###  II.2 Wave amplitude equation

Now, we introduce the normal variables $a^{s}(\mathbf{k},t)$ defined as
follows:

$a^{s}(\mathbf{k},t)=\frac{1}{\epsilon}\left(\sqrt{\frac{k}{2}}\hat{\varphi}(\mathbf{k},t)+\frac{is}{\sqrt{2k}}\partial_{t}\hat{\varphi}(\mathbf{k},t)\right)e^{is\omega_{\mathbf{k}}t},$
(8)

with $s=\pm 1$ the directional polarity (for a given $\mathbf{k}$ we have two
directions of propagation) and $\hat{\varphi}(\mathbf{k},t)$ the spatial
Fourier transform of $\varphi(\mathbf{x},t)$. Conversely, we have the
relations:

$\hat{\varphi}(\mathbf{k},t)=\frac{\epsilon}{\sqrt{2k}}\sum_{s}a^{s}_{k}e^{-is\omega_{\mathbf{k}}t},$
(9a) and
$\partial_{t}\hat{\varphi}(\mathbf{k},t)=-i\epsilon\sqrt{\frac{k}{2}}\sum_{s}sa^{s}_{k}e^{-is\omega_{\mathbf{k}}t}.$
(9b)

Referring to the supplementary material of [20], we obtain explicit
expressions for $\alpha(\mathbf{k},t)$, $\beta(\mathbf{k},t)$, and
$\gamma(\mathbf{k},t)$ in terms of the normal variables using equations (7a),
(7b), and (7c). In this context, we assume that the amplitude evolves much
slower than the phase, namely:

$\frac{\partial_{t}a^{s}(\mathbf{k})}{a^{s}(\mathbf{k})}\ll\omega_{\bm{k}}.$
(10)

This assumption, which is fully consistent with the multiple timescale method
introduced in the next section, allows us to treat the amplitude as
approximately constant compared to the phase variation during a time
integration. Additionally, we consider the integration variables as “dumb”,
enabling us to symmetrize the integrands.

The Fourier transform of equations (7a)–(7b) leads to the following
expressions:

$\dot{\alpha}(\mathbf{k})=i\frac{\epsilon^{2}}{2}\int_{\mathbb{R}^{4}}\sum_{s_{1},s_{2}}\frac{s_{1}k_{1}p_{2}+s_{2}k_{2}p_{1}}{p\sqrt{k_{1}k_{2}}}a^{s_{1}}(\mathbf{k}_{1})a^{s_{2}}(\mathbf{k}_{2})e^{i\Omega_{12}t}\delta^{0}_{12}(\mathbf{k})\mathrm{d}^{2}\mathbf{k}_{1}\mathrm{d}^{2}\mathbf{k}_{2},$
(11)

and

$\dot{\beta}(\mathbf{k})=i\frac{\epsilon^{2}}{2}\int_{\mathbb{R}^{4}}\sum_{s_{1},s_{2}}\frac{s_{1}k_{1}q_{2}+s_{2}k_{2}q_{1}}{q\sqrt{k_{1}k_{2}}}a^{s_{1}}(\mathbf{k}_{1})a^{s_{2}}(\mathbf{k}_{2})e^{i\Omega_{12}t}\delta^{0}_{12}(\mathbf{k})\mathrm{d}^{2}\mathbf{k}_{1}\mathrm{d}^{2}\mathbf{k}_{2}.$
(12)

These equations can be integrated over time, resulting in the following
expressions:

$\alpha(\mathbf{k})=-\frac{\epsilon^{2}}{2}\int_{\mathbb{R}^{4}}\sum_{s_{1},s_{2}}\frac{s_{1}k_{1}p_{2}+s_{2}k_{2}p_{1}}{p\sqrt{k_{1}k_{2}}}\frac{a^{s_{1}}(\mathbf{k}_{1})a^{s_{2}}(\mathbf{k}_{2})}{s_{1}k_{1}+s_{2}k_{2}}e^{i\Omega_{12}t}\delta^{0}_{12}(\mathbf{k})\mathrm{d}^{2}\mathbf{k}_{1}\mathrm{d}^{2}\mathbf{k}_{2}$
(13)

and,

$\beta(\mathbf{k})=-\frac{\epsilon^{2}}{2}\int_{\mathbb{R}^{4}}\sum_{s_{1},s_{2}}\frac{s_{1}k_{1}q_{2}+s_{2}k_{2}q_{1}}{q\sqrt{k_{1}k_{2}}}\frac{a^{s_{1}}(\mathbf{k}_{1})a^{s_{2}}(\mathbf{k}_{2})}{s_{1}k_{1}+s_{2}k_{2}}e^{i\Omega_{12}t}\delta^{0}_{12}(\mathbf{k})\mathrm{d}^{2}\mathbf{k}_{1}\mathrm{d}^{2}\mathbf{k}_{2}.$
(14)

In addition, the Fourier transform of equation (7c) leads to the following
expression:

$\gamma(\mathbf{k})=-\frac{\epsilon^{2}}{2}\int_{\mathbb{R}^{4}}\sum_{s_{1},s_{2}}\frac{p_{1}q_{2}+p_{2}q_{1}}{pq}\frac{1}{\sqrt{k_{1}k_{2}}}a^{s_{1}}(\mathbf{k}_{1})a^{s_{2}}(\mathbf{k}_{2})e^{i\Omega_{12}t}\delta^{0}_{12}(\mathbf{k})\mathrm{d}^{2}\mathbf{k}_{1}\mathrm{d}^{2}\mathbf{k}_{2}.$
(15)

This expression can be differentiated with respect to time, resulting in:

$\dot{\gamma}(\mathbf{k})=i\frac{\epsilon^{2}}{2}\int_{\mathbb{R}^{4}}\sum_{s_{1},s_{2}}\frac{p_{1}q_{2}+p_{2}q_{1}}{pq}\frac{s_{1}k_{1}+s_{2}k_{2}}{\sqrt{k_{1}k_{2}}}a^{s_{1}}(\mathbf{k}_{1})a^{s_{2}}(\mathbf{k}_{2})e^{i\Omega_{12}t}\delta^{0}_{12}(\mathbf{k})\mathrm{d}^{2}\mathbf{k}_{1}\mathrm{d}^{2}\mathbf{k}_{2}.$
(16)

Next, we substitute all the previously derived expressions into equation (7d)
in order to find the wave amplitude equation. We obtain after manipulations:

$\partial_{t}a^{s}(\mathbf{k})=\epsilon^{2}\int_{\mathbb{R}^{6}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}a^{s_{1}}(\mathbf{k}_{1})a^{s_{2}}(\mathbf{k}_{2})a^{s_{3}}(\mathbf{k}_{3})e^{i\Omega^{0}_{123}t}~{}\delta^{0}_{123}(\mathbf{k})~{}\prod_{i=1}^{3}\mathrm{d}^{2}\mathbf{k}_{i},$
(17)

where
$\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}$
is the interaction coefficient, defined as:

$\mathbf{L}^{ss_{1}s_{2}s_{3}}_{\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}=\frac{is}{12\sqrt{kk_{1}k_{2}k_{3}}}\left(\mathbf{T}^{ss_{1}s_{2}s_{3}}_{\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}+\mathbf{T}^{ss_{2}s_{3}s_{1}}_{\mathbf{k}\mathbf{k}_{2}\mathbf{k}_{3}\mathbf{k}_{1}}+\mathbf{T}^{ss_{3}s_{2}s_{1}}_{\mathbf{k}\mathbf{k}_{3}\mathbf{k}_{2}\mathbf{k}_{1}}\right),$
(18)

with:

$\displaystyle\mathbf{T}^{ss_{1}s_{2}s_{3}}_{\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}=$
$\displaystyle\bigg{[}\left(p_{2}q_{3}+p_{3}q_{2}\right)\frac{pp_{1}+qq_{1}-s_{1}k_{1}\left(s_{1}k_{1}+s_{2}k_{2}+s_{3}k_{3}\right)}{\left(p_{2}+p_{3}\right)\left(q_{2}+q_{3}\right)}$
(19) $\displaystyle\hskip
5.0pt+\frac{s_{2}k_{2}p_{3}+s_{3}k_{3}p_{2}}{s_{2}k_{2}+s_{3}k_{3}}\frac{pp_{1}-qq_{1}+s_{1}k_{1}\left(s_{1}k_{1}+s_{2}k_{2}+s_{3}k_{3}\right)}{p_{2}+p_{3}}$
$\displaystyle\hskip
5.0pt-\frac{s_{2}k_{2}q_{3}+s_{3}k_{3}q_{2}}{s_{2}k_{2}+s_{3}k_{3}}\frac{pp_{1}-qq_{1}-s_{1}k_{1}\left(s_{1}k_{1}+s_{2}k_{2}+s_{3}k_{3}\right)}{q_{2}+q_{3}}\bigg{]}.$

###  II.3 Symmetries of the interaction coefficient

The specific form of the interaction coefficient is determined by considering
the dummy character of the integration variable, which allows us to identify
certain symmetries of the coefficient. In general, the interaction coefficient
exhibits the following properties:

  *   * •

Homogeneity: It is a homogeneous function of the wave vectors
$\mathbf{k}_{i}$. This means that there exists a constant $b\in\mathbb{R}$
such that, for all $a\in\mathbb{R}$,

$\mathbf{L}^{ss_{1}s_{2}s_{3}}_{a\mathbf{k}a\mathbf{k}_{1}a\mathbf{k}_{2}a\mathbf{k}_{3}}=a^{b}\mathbf{L}^{ss_{1}s_{2}s_{3}}_{\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}.$
(20a)

It is worth noting that, for our specific problem, $b$ is equal to zero.

  * •

Imaginary nature: The interaction coefficient is purely imaginary, satisfying
the relation:

$\left(\mathbf{L}^{ss_{1}s_{2}s_{3}}_{\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\right)^{*}=-\mathbf{L}^{ss_{1}s_{2}s_{3}}_{\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}.$
(20b)

  * •

Symmetry: The interaction coefficient exhibits symmetry under certain
transformations. Specifically, it is symmetric when the signs of the
directional polarities or all wave vectors are changed:

$\mathbf{L}^{-s-s_{1}-s_{2}-s_{3}}_{-\mathbf{k}-\mathbf{k}_{1}-\mathbf{k}_{2}-\mathbf{k}_{3}}=-\mathbf{L}^{ss_{1}s_{2}s_{3}}_{\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}.$
(20c)

  * •

Symmetries on the resonant manifold: On the resonant manifold, where
$s\omega_{\mathbf{k}}=s_{1}\omega_{\mathbf{k}_{1}}+s_{2}\omega_{\mathbf{k}_{2}}+s_{3}\omega_{\mathbf{k}_{3}}$
and $\mathbf{k}=\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}$, the interaction
coefficient exhibits symmetry with respect to the last indices. Thus, for any
permutation $\sigma$ of $\left\\{1,2,3\right\\}$, we have:

$\mathbf{L}^{ss_{\sigma(1)}s_{\sigma(2)}s_{\sigma(3)}}_{\mathbf{k}\mathbf{k}{\sigma(1)}\mathbf{k}{\sigma(2)}\mathbf{k}{\sigma(3)}}=\mathbf{L}^{ss_{1}s_{2}s_{3}}_{\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}.$
(20d)

Additionally, the interaction coefficient vanishes when $\mathbf{k}=0$,
ensuring the conservation of the null first-order cumulant:

$\mathbf{L}^{ss_{1}s_{2}s_{3}}_{\mathbf{0}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}=0.$
(20e)

##  III Sequential time closure for quartic-wave interactions

###  III.1 A multiple time scale method

Our statistical study begins with the fundamental equation (17) from which we
apply the method of multiple time scales to derive the kinetic equation. Note
that this method was originally developed for three-wave interactions [22] and
then to quartic interactions with gravity waves as its main application [24,
1]. The original problem (surface waves) is 2-dimensional, however, the method
is valid in a space of n-dimension, with $n\geq 1$. That is why, for this
section only, we will perform the spatial integrations on vectors of
$\mathbb{R}^{n}$.

We introduce a set of (approximately) independent variables:

$T_{0}=t,\quad T_{1}=\epsilon t,\quad T_{2}=\epsilon^{2}t,\quad\dots$ (21)

with $\epsilon$ a small parameter ($0<\epsilon\ll 1$) which measures the wave
amplitude (see section II). We apply the chain rule, yielding:

$\partial_{t}=\sum_{i=0}^{+\infty}\epsilon^{i}\partial_{T_{i}}=\partial_{T_{0}}+\epsilon\partial_{T_{1}}+\epsilon^{2}\partial_{T_{2}}+\dots$
(22)

The normal variables are also expanded to the power of $\epsilon$ as follows:

$a^{s}(\mathbf{k})=\sum_{i=0}^{+\infty}\epsilon^{i}a^{s}_{i}(\mathbf{k})=a^{s}_{0}(\mathbf{k})+\epsilon
a^{s}_{1}(\mathbf{k})+\epsilon^{2}a^{s}_{2}(\mathbf{k})+\dots$ (23)

In order to enhance readability, we disregard the time dependence, so
$a^{s}(\mathbf{k})$ must be understood like
$a^{s}(\mathbf{k},T_{0},T_{1},T_{2},\dots)$. Furthermore, we assume that the
phase of the normal variables only varies on the $T_{0}$ timescale. The
justification for this assumption will be provided later.

With these definitions, we obtain the following expansion:

$\sum_{i,j=0}^{+\infty}\epsilon^{i+j}\partial_{T_{i}}a^{s}_{j}(\mathbf{k})=\\\
\sum_{p,q,r=0}^{+\infty}\epsilon^{p+q+r+2}\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}a^{s_{1}}_{p}(\mathbf{k}_{1})a^{s_{2}}_{q}(\mathbf{k}_{2})a^{s_{3}}_{r}(\mathbf{k}_{3})~{}\delta^{0}_{123}(\mathbf{k})~{}e^{i\Omega^{0}_{123}T_{0}}~{}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}.$
(24)

This expression leads to a hierarchy of equations in different orders in
$\epsilon$. Using eq. (24) to identify the different powers of $\epsilon$, we
obtain the following expressions:

$\displaystyle\partial_{T_{0}}a^{s}_{0}(\mathbf{k})=$ $\displaystyle~{}0;$
(25a) $\displaystyle\partial_{T_{0}}a^{s}_{1}(\mathbf{k})=$
$\displaystyle-\partial_{T_{1}}a^{s}_{0}(\mathbf{k});$ (25b)
$\displaystyle\partial_{T_{0}}a^{s}_{2}(\mathbf{k})=$
$\displaystyle-\partial_{T_{2}}a^{s}_{0}(\mathbf{k})-\partial_{T_{1}}a^{s}_{1}(\mathbf{k})$
(25c)
$\displaystyle+\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}a^{s_{1}}_{0}(\mathbf{k}_{1})a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})\delta^{0}_{123}(\mathbf{k})e^{i\Omega^{0}_{123}T_{0}}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i};$
(25d) $\displaystyle\partial_{T_{0}}a^{s}_{3}(\mathbf{k})=$
$\displaystyle-\partial_{T_{3}}a^{s}_{0}(\mathbf{k})-\partial_{T_{2}}a^{s}_{1}(\mathbf{k})-\partial_{T_{1}}a^{s}_{2}(\mathbf{k})$
(25e)
$\displaystyle+3\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}a^{s_{1}}_{1}(\mathbf{k}_{1})a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})\delta^{0}_{123}(\mathbf{k})e^{i\Omega^{0}_{123}T_{0}}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i};$
(25f)
$\displaystyle\begin{split}\partial_{T_{0}}a^{s}_{4}(\mathbf{k})=&-\partial_{T_{4}}a^{s}_{0}(\mathbf{k})-\partial_{T_{3}}a^{s}_{1}(\mathbf{k})-\partial_{T_{2}}a^{s}_{2}(\mathbf{k})-\partial_{T_{1}}a^{s}_{3}(\mathbf{k})\\\
&+3\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}a^{s_{1}}_{2}(\mathbf{k}_{1})a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})\delta^{0}_{123}(\mathbf{k})e^{i\Omega^{0}_{123}T_{0}}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}\\\
&+3\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}a^{s_{1}}_{1}(\mathbf{k}_{1})a^{s_{2}}_{1}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})\delta^{0}_{123}(\mathbf{k})e^{i\Omega^{0}_{123}T_{0}}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}.\end{split}$
(25g)

The factor of $3$ arises from considering the symmetries of the interaction
coefficient. Additionally, we observe that eq. (25a) provides justification
for assuming that the amplitude and phase do not vary on the same timescale.

We choose initial conditions (at $t=0$) such that: $\forall i\geq 1$,
$a_{i}^{s}(\mathbf{k})=0$. The integration of the previous equation with
respect to $T_{0}$ gives, after a few manipulations:

$\displaystyle a^{s}_{0}(\mathbf{k})=$
$\displaystyle~{}a^{s}_{0}(\mathbf{k},T_{1},T_{2},\dots);$ (26a)
$\displaystyle a^{s}_{1}(\mathbf{k})=$ $\displaystyle-
T_{0}\partial_{T_{1}}a^{s}_{0}(\mathbf{k});$ (26b) $\displaystyle
a^{s}_{2}(\mathbf{k})=$ $\displaystyle-
T_{0}\partial_{T_{2}}a^{s}_{0}(\mathbf{k})+\frac{T_{0}^{\phantom{a}2}}{2}\partial^{2}_{T_{1}}a^{s}_{0}(\mathbf{k})+b^{s}_{2}(\mathbf{k});$
(26c) $\displaystyle a^{s}_{3}(\mathbf{k})=$ $\displaystyle-
T_{0}\partial_{T_{3}}a^{s}_{0}(\mathbf{k})+T_{0}^{\phantom{a}2}\partial_{T_{1}}\partial_{T_{2}}a^{s}_{0}(\mathbf{k})-\frac{T_{0}^{\phantom{a}3}}{6}\partial^{3}_{T_{1}}a^{s}_{0}(\mathbf{k})-T_{0}\partial_{T_{1}}b_{2}^{s}(\mathbf{k});$
(26d)
$\displaystyle\begin{split}a^{s}_{4}(\mathbf{k})=&-T_{0}\partial_{T_{4}}a^{s}_{0}(\mathbf{k})+\frac{T_{0}^{\phantom{a}2}}{2}\left[2\partial_{T_{1}}\partial_{T_{3}}a^{s}_{0}(\mathbf{k})+\partial_{T_{2}}^{2}a^{s}_{0}(\mathbf{k})\right]-\frac{T_{0}^{\phantom{a}3}}{2}\partial_{T_{1}}^{2}\partial_{T_{2}}a^{s}_{0}(\mathbf{k})\\\
&+\frac{T_{0}^{\phantom{a}4}}{24}\partial_{T_{1}}^{4}a^{s}_{0}(\mathbf{k})+b_{4}^{s}(\mathbf{k})-T_{0}\partial_{T_{2}}b_{2}^{s}(\mathbf{k})+\frac{T_{0}^{\phantom{a}2}}{2}\partial_{T_{1}}^{2}b_{2}^{s}(\mathbf{k});\end{split}$
(26e)

where:

$\displaystyle\displaystyle b^{s}_{2}(\mathbf{k})=$
$\displaystyle~{}\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}a^{s_{1}}_{0}(\mathbf{k}_{1})a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})\delta^{0}_{123}(\mathbf{k})\Delta\left(\Omega^{0}_{123},T_{0}\right)\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i};$
(27a)
$\displaystyle\begin{split}b^{s}_{4}(\mathbf{k})=&~{}3\int_{\mathbb{R}^{6n}}\sum_{\begin{subarray}{c}s_{1},s_{2},s_{3}\\\
s_{4},s_{5},s_{6}\end{subarray}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\mathbf{L}^{-s_{1}s_{4}s_{5}s_{6}}_{-\mathbf{k}_{1}\mathbf{k}_{4}\mathbf{k}_{5}\mathbf{k}_{6}}a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})a^{s_{4}}_{0}(\mathbf{k}_{4})a^{s_{5}}_{0}(\mathbf{k}_{5})a^{s_{6}}_{0}(\mathbf{k}_{6})\\\\[-18.0pt]
&\hskip
131.0pt\times\delta^{0}_{123}(\mathbf{k})\delta^{1}_{456}(\mathbf{k})E\left(\Omega^{0}_{23456},\Omega^{0}_{123},T_{0}\right)\prod_{i=1}^{6}\mathrm{d}^{n}\mathbf{k}_{i};\end{split}$
(27b)

with:

$\Delta\left(X,T\right)=\frac{e^{iXT}-1}{iX}\quad\text{and}\quad
E\left(X,Y,T\right)=\int_{0}^{T}\Delta(X-Y,t)e^{iYt}~{}\mathrm{d}t.$

###  III.2 Statistical assumptions

We make the statistical assumption that $\langle a^{s}(\mathbf{k})\rangle=0$,
where $\langle\rangle$ is the ensemble average. Going back to general
relativity, this situation corresponds to space-time fluctuations around a
flat metric. We also assume that $\varphi(\mathbf{r})$ and
$\partial_{t}\varphi(\mathbf{r})$ are spatially homogeneous, meaning that we
are in the heart of the turbulence. As a result, the corresponding normal
variables $a^{s}(\mathbf{k})$ are also spatially homogeneous since they are
defined as linear combinations of these fields. The second-order moment can
thus be written in terms of the second-order cumulant
$q^{ss^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})$ as follows:

$\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle=q_{0}^{ss^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})\delta\left(\mathbf{k}+\mathbf{k}^{\prime}\right),$
(28)

where the presence of $\delta\left(\mathbf{k}+\mathbf{k}^{\prime}\right)$ is
the consequence of statistical homogeneity. Note that for the derivation of
the kinetic equation, the product $\langle
a^{s}(\mathbf{k})a^{s^{\prime}}(\mathbf{k}^{\prime})\rangle$ will finally be
treated in the special case where $s^{\prime}=-s$. The reason is that the
quantities of interest are wave action
($\langle|a^{s}(\mathbf{k})|^{2}\rangle=\langle
a^{s}(\mathbf{k})a^{-s}(\mathbf{-k})\rangle$) and energy
($\omega_{\mathbf{k}}\langle\left|a^{s}(\mathbf{k})\right|^{2}\rangle$)
spectra which are real and positive.

The $n$th-order moment $\langle a^{s_{1}}(\mathbf{k}_{1})\dots
a^{s_{r}}(\mathbf{k}_{r})\rangle$ can be written in terms of lower order
cumulants $q^{s_{1}\dots s_{r}}(\mathbf{k}_{1},\dots,\mathbf{k}_{r})$. For
instance, if $r=4$, we have for the leading order cumulants:

$\displaystyle\langle
a^{s_{1}}_{0}(\mathbf{k}_{1})a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})a^{s_{4}}_{0}(\mathbf{k}_{4})\rangle=$
$\displaystyle~{}q_{0}^{s_{1}s_{2}s_{3}s_{4}}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4})\delta^{1234}(\mathbf{k})$
(29)
$\displaystyle+q_{0}^{s_{1}s_{2}}(\mathbf{k}_{1},\mathbf{k}_{2})q_{0}^{s_{3}s_{4}}(\mathbf{k}_{3},\mathbf{k}_{4})\delta^{12}(\mathbf{k})\delta^{34}(\mathbf{k})$
$\displaystyle+q_{0}^{s_{1}s_{3}}(\mathbf{k}_{1},\mathbf{k}_{3})q_{0}^{s_{2}s_{4}}(\mathbf{k}_{2},\mathbf{k}_{4})\delta^{13}(\mathbf{k})\delta^{24}(\mathbf{k})$
$\displaystyle+q_{0}^{s_{1}s_{4}}(\mathbf{k}_{1},\mathbf{k}_{4})q_{0}^{s_{2}s_{3}}(\mathbf{k}_{2},\mathbf{k}_{3})\delta^{14}(\mathbf{k})\delta^{23}(\mathbf{k}).$

Note that the useful case $r=6$ is recalled in appendix A.

To simplify the analysis, we will assume Gaussian 111Note that the initial
Gaussian assumption is a priori not required for the derivation of the kinetic
equation, and we can also assume, as often [33], initial fields
$a_{0}^{s}(\mathbf{k})$ with random phase and amplitude. initial distribution
for $a_{0}^{s}(\mathbf{k})$ so that the fourth-order cumulant at order zero
vanishes,
$q^{ss^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}_{0}(\mathbf{k},\mathbf{k}^{\prime},\mathbf{k}^{\prime\prime},\mathbf{k}^{\prime\prime\prime},t=0)=0$.
This assumption is sufficient to guarantee the absence of coherent structure
in the initial excitation that could invalidate our derivation based on weak
nonlinearities. Under the additional assumption of boundedness of the fourth-
order moment when $T_{0}\to+\infty$, we will show that this initial property
propagates at time $T_{2}$ and even has an impact on the kinetic equation at
time $T_{4}$. This phenomenon is sometimes referred to as propagation of chaos
[44].

A fundamental aspect of the multiple timescale method is the assumption that
the second-order moment $\langle
a^{s}(\mathbf{k})a^{s^{\prime}}(\mathbf{k}^{\prime})\rangle$ remains bounded
as $T_{0}\to+\infty$. Making use of the previous definitions, we expend this
term according to:

$\displaystyle\langle
a^{s}(\mathbf{k})a^{s^{\prime}}(\mathbf{k}^{\prime})\rangle$
$\displaystyle=\left\langle\left(a^{s}_{0}(\mathbf{k})+\epsilon
a^{s}_{1}(\mathbf{k})+\epsilon^{2}a^{s\phantom{{}^{\prime}}}_{2}(\mathbf{k})+\dots\right)\right.$
(30)
$\displaystyle\times\left.\left(a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})+\epsilon
a^{s^{\prime}}_{1}(\mathbf{k}^{\prime})+\epsilon^{2}a^{s^{\prime}}_{2}(\mathbf{k}^{\prime})+\dots\right)\right\rangle$
$\displaystyle=\sum_{n=0}^{+\infty}\epsilon^{n}\sum_{i=0}^{n}\langle
a^{s}_{i}(\mathbf{k})a^{s}_{n-i}(\mathbf{k}^{\prime})\rangle.$

Thus, the contribution of the right-hand side must also be bounded at each
order $\epsilon^{n}$, with $n\geq 0$. In practice, we will see that secular
terms may appear and, therefore, we will impose the nullity of these
contributions at a given order to keep the development uniform in time. It is
precisely this condition at order $\epsilon^{4}$ that will lead us to the
kinetic equation.

Our last statistical assumption concerns the boundedness of the fourth-order
cumulants over time. In other words, we impose that $\langle
a^{s_{1}}(\mathbf{k}_{1})a^{s_{2}}(\mathbf{k}_{2})a^{s_{3}}(\mathbf{k}_{3})a^{s_{4}}(\mathbf{k}_{4})\rangle$
does not diverge when $T_{0}\to+\infty$. Thus, we will have to consider sums
of products given by:

$\sum_{\begin{subarray}{c}i_{1},i_{2},i_{3},i_{4}\\\
i_{1}+i_{2}+i_{3}+i_{4}=n\end{subarray}}\langle
a^{s_{1}}_{i_{1}}(\mathbf{k}_{1})a^{s_{2}}_{i_{2}}(\mathbf{k}_{2})a^{s_{3}}_{i_{3}}(\mathbf{k}_{3})a^{s_{4}}_{i_{4}}(\mathbf{k}_{4})\rangle,$
(31)

and impose the nullity of the corresponding secular terms.

###  III.3 First asymptotic closure at time $\mathbf{T_{1}}$

Now we will derive the conditions which ensure the second-order moment remains
bounded at all times. This involves counting all the secular terms at each
order $\epsilon^{n}$, $n\geq 0$ and cancelling their contributions. We will
see that it is the condition of order $\epsilon^{4}$ that gives the kinetic
equation we are looking for. Using expression (30), we see that at order
$\epsilon^{0}$ the result is immediate: the only term is $\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle$ which is
bounded, so there is no secular contribution and therefore no asymptotic
condition to impose.

At order $\epsilon$ the situation is different because we have:

$\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{1}(\mathbf{k}^{\prime})+a^{s}_{1}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle=-T_{0}\partial_{T_{1}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle.$ (32)

In order to keep the left-hand side bounded, we must impose the condition:

$\partial_{T_{1}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle=0,$ (33)

which means that the wave action does not evolve on the $T_{1}$ timescale.
Therefore, the turbulent cascade is not effective on this short time scale.
Similarly, we can justify that:

$\partial_{T_{1}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})a^{s^{\prime\prime}}_{0}(\mathbf{k}^{\prime\prime})a^{s^{\prime\prime\prime}}_{0}(\mathbf{k}^{\prime\prime\prime})\rangle=0,$
(34)

based on the bounded nature of the fourth-order moment222These results based
on the boundedness of second and fourth order moments can be assumed for any
moments of order $n$. Therefore, the probability distribution does not depend
on $T_{1}$ and we may assume that the $a^{s}_{0}(\mathbf{k})$ variable itself
does not depend on $T_{1}$ either. Thus, at most we have:
$a^{s}_{0}(\mathbf{k})=a^{s}_{0}(\mathbf{k},T_{2},T_{3},\dots)$. This
justifies the requirement that the $T_{1}$ timescale is irrelevant in the
derivation of the kinetic equation as it was done by Newell [1]. Nevertheless,
in our derivation, we will keep the $T_{1}$ dependence to show that this
dependence can disappear with this boundedness condition.

###  III.4 Second asymptotic closure at time $\mathbf{T_{2}}$

We continue the analysis with the second order:

$\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{2}(\mathbf{k}^{\prime})+a^{s}_{1}(\mathbf{k})a^{s^{\prime}}_{1}(\mathbf{k}^{\prime})+a^{s}_{2}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle\\\
=-T_{0}\partial_{T_{2}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle+\frac{T_{0}^{\phantom{a}2}}{2}\partial_{T_{1}}^{2}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle+\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})+a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})b^{s}_{2}(\mathbf{k})\rangle.$
(35)

This expression can be simplified using the first closure (33). Therefore, to
prevent unwanted secular growth of the left-hand side, we must balance the
terms proportional to $T_{0}$:

$\partial_{T_{2}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle=\mathcal{C}_{T_{0}}\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})+a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})b^{s}_{2}(\mathbf{k})\rangle,$
(36)

where $\mathcal{C}_{T_{0}^{\phantom{a}n}}\mathcal{E}$ refers to the term
proportional to $T_{0}^{\phantom{a}n}$ in $\mathcal{E}$. The long time
contribution of these oscillating integrals is given by the theory of
generalized function and the Riemann-Lebesgue lemma (see Appendix B):

$\Delta(X,T)\xrightarrow[T_{0}\to+\infty]{}\begin{cases}\pi\delta(X)+i\mathcal{P}\left(\frac{1}{X}\right),&\text{if}~{}X\neq
0\\\ T_{0},&\text{if}~{}X=0\\\ \end{cases};$ (37)

where $\mathcal{P}$ refers to the Cauchy principal value. This means that
$\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})+a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})b^{s}_{2}(\mathbf{k})\rangle$
can exhibit secular growth when the different Dirac deltas constrain the
integration to be performed on the resonant manifold.

We illustrate the computation with an example:

$\displaystyle\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})\rangle$ (38)
$\displaystyle=\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-s^{\prime}s_{1}s_{2}s_{3}}_{-\mathbf{k}^{\prime}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\langle
a^{s}_{0}(\mathbf{k})a^{s_{1}}_{0}(\mathbf{k}_{1})a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})\rangle\Delta(\Omega^{0^{\prime}}_{123},T_{0})\delta^{0^{\prime}}_{123}(\mathbf{k})\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}$
$\displaystyle=~{}3\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-s^{\prime}s_{1}s_{2}s_{3}}_{-\mathbf{k}^{\prime}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}q_{0}^{ss_{1}}(\mathbf{k},\mathbf{k}_{1})q_{0}^{s_{2}s_{3}}(\mathbf{k}_{2},\mathbf{k}_{3})\delta^{01}(\mathbf{k})\delta^{23}(\mathbf{k})\Delta(\Omega^{0^{\prime}}_{123},T_{0})\delta^{0^{\prime}}_{123}(\mathbf{k})\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}$
$\displaystyle+\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-s^{\prime}s_{1}s_{2}s_{3}}_{-\mathbf{k}^{\prime}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}q_{0}^{ss_{1}s_{2}s_{3}}(\mathbf{k},\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})\delta^{0123}(\mathbf{k})\Delta(\Omega^{0^{\prime}}_{123},T_{0})\delta^{0^{\prime}}_{123}(\mathbf{k})\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}.$

The first term can exhibit secular growth if we select the polarization as
$s_{1}=s^{\prime}$ and $s_{2}=-s_{3}$. In this case, the Dirac delta functions
ensure that $\omega_{\mathbf{k}^{\prime}}=\omega_{\mathbf{k}_{1}}$ and
$\omega_{\mathbf{k}_{2}}=\omega_{\mathbf{k}_{3}}$, which leads to the
vanishing of $\Omega^{0^{\prime}}_{123}$. The remaining terms do not
contribute to secular growth. Therefore, we obtain the following result:

$\mathcal{C}_{T_{0}}\left\\{\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})\rangle\right\\}=3\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle
F(\mathbf{k}^{\prime},s^{\prime}),$ (39)

where:

$F(\mathbf{k}^{\prime},s^{\prime})=\int_{\mathbb{R}^{n}}\sum_{s_{2}}\mathbf{L}^{-s^{\prime}s^{\prime}s_{2}-s_{2}}_{-\mathbf{k}^{\prime}\mathbf{k}^{\prime}\mathbf{k}_{2}-\mathbf{k}_{2}}q_{0}^{s_{2}-s_{2}}(\mathbf{k}_{2},-\mathbf{k}_{2})~{}\mathrm{d}^{n}\mathbf{k}_{2}.$
(40)

Combining both secular contributions, we have:

$\partial_{T_{2}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle=\\\
3\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle\left[F(\mathbf{k},s)+F(\mathbf{k}^{\prime},s^{\prime})\right].$
(41)

However, recalling that we are in the case where
$\mathbf{k}^{\prime}=-\mathbf{k}$ and $s^{\prime}=-s$ which is the only
relevant case for our study, we may notice that:

$\displaystyle F(\mathbf{k}^{\prime},s^{\prime})$
$\displaystyle=\int_{\mathbb{R}^{n}}\sum_{s_{2}}\mathbf{L}^{s-ss_{2}-s_{2}}_{\mathbf{k}-\mathbf{k}\mathbf{k}_{2}-\mathbf{k}_{2}}q_{0}^{s_{2}-s_{2}}(\mathbf{k}_{2},-\mathbf{k}_{2})~{}\mathrm{d}^{n}\mathbf{k}_{2}$
(42)
$\displaystyle=-\int_{\mathbb{R}^{n}}\sum_{s_{2}}\mathbf{L}^{-sss_{2}-s_{2}}_{-\mathbf{k}\mathbf{k}\mathbf{k}_{2}-\mathbf{k}_{2}}q_{0}^{s_{2}-s_{2}}(\mathbf{k}_{2},-\mathbf{k}_{2})~{}\mathrm{d}^{n}\mathbf{k}_{2}$
$\displaystyle=-F(\mathbf{k},s)$

which implies that $\partial_{T_{2}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle$ vanishes
when we consider the symmetries of the interaction, resulting in:

$\partial_{T_{2}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle=0.$ (43)

Once again, the wave action does not evolve with respect to $T_{2}$ and the
turbulent cascade does not develop at this short timescale.

This analysis can be extended to the fourth-order moment. We obtain:

$\partial_{T_{2}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})a^{s^{\prime\prime}}_{0}(\mathbf{k}^{\prime\prime})a^{s^{\prime\prime\prime}}_{0}(\mathbf{k}^{\prime\prime\prime})\rangle\\\
=3\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})a^{s^{\prime\prime}}_{0}(\mathbf{k}^{\prime\prime})a^{s^{\prime\prime\prime}}_{0}(\mathbf{k}^{\prime\prime\prime})\rangle\left[F(\mathbf{k},s)+F(\mathbf{k}^{\prime},s^{\prime})+F(\mathbf{k}^{\prime\prime},s^{\prime\prime})+F(\mathbf{k}^{\prime\prime\prime},s^{\prime\prime\prime})\right].$
(44)

From the second closure, we know that the second-order moment does not
depend333This result has been proved when $s^{\prime}=-s$. In the other case,
where $s^{\prime}=s$ the statistical mean $\langle
a^{s}(\mathbf{k},t)a^{s^{\prime}}(\mathbf{k}^{\prime},t)e^{i(s\omega_{\mathbf{k}}+s^{\prime}\omega_{\mathbf{k}^{\prime}})t}\rangle=q(\mathbf{k},\mathbf{k}^{\prime})\delta(\mathbf{k}+\mathbf{k}^{\prime})\langle
e^{i(s\omega_{\mathbf{k}}+s^{\prime}\omega_{\mathbf{k}^{\prime}})t}\rangle$ is
necessary zero since the dirac function impose that
$\omega_{\mathbf{k}}=\omega_{\mathbf{k}^{\prime}}$ so oscillating term
vanishes. Thus, the case $s^{\prime}=s$ has no relevant physical meaning. on
$T_{2}$, and therefore the left hand side can be reduced to a fourth-order
cumulant. On the right hand side, we see that the decomposition in terms of
second-order moments does not give any contribution because of the pairwise
cancellation. Therefore, the previous equation can be reduced to:

$\partial_{T_{2}}q_{0}^{ss^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(\mathbf{k},\mathbf{k^{\prime}},\mathbf{k^{\prime\prime}},\mathbf{k^{\prime\prime\prime}})\\\
=3q_{0}^{ss^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(\mathbf{k},\mathbf{k^{\prime}},\mathbf{k^{\prime\prime}},\mathbf{k^{\prime\prime\prime}})\left[F(\mathbf{k},s)+F(\mathbf{k}^{\prime},s^{\prime})+F(\mathbf{k}^{\prime\prime},s^{\prime\prime})+F(\mathbf{k}^{\prime\prime\prime},s^{\prime\prime\prime})\right].$
(45)

This relation can be integrated according to $T_{2}$ timescale, giving:

$q_{0}^{ss^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(\mathbf{k},\mathbf{k^{\prime}},\mathbf{k^{\prime\prime}},\mathbf{k^{\prime\prime\prime}})\\\
=\left.q_{0}^{ss^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(\mathbf{k},\mathbf{k^{\prime}},\mathbf{k^{\prime\prime}},\mathbf{k^{\prime\prime\prime}})\right|_{t=0}\exp\left\\{3T_{2}\left[F(\mathbf{k},s)+F(\mathbf{k}^{\prime},s^{\prime})+F(\mathbf{k}^{\prime\prime},s^{\prime\prime})+F(\mathbf{k}^{\prime\prime\prime},s^{\prime\prime\prime})\right]\right\\}.$
(46)

However, according to the initial conditions (Gaussian statistics or random
phase assumption)

$\left.q_{0}^{ss^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(\mathbf{k},\mathbf{k^{\prime}},\mathbf{k^{\prime\prime}},\mathbf{k^{\prime\prime\prime}})\right|_{t=0}=0,$
(47)

which implies that:

$q_{0}^{ss^{\prime}s^{\prime\prime}s^{\prime\prime\prime}}(\mathbf{k},\mathbf{k^{\prime}},\mathbf{k^{\prime\prime}},\mathbf{k^{\prime\prime\prime}})=0,$
(48)

and its first derivative in $T_{2}$ is also equal to zero. Therefore, we have
proved that this Gaussian/random phase property assumed initially propagates
(at least) up to $T_{2}$ when initially assumed.

###  III.5 Third asymptotic closure at time $\mathbf{T_{3}}$

At next order, we have:

$\displaystyle\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{3}(\mathbf{k}^{\prime})+a^{s}_{1}(\mathbf{k})a^{s^{\prime}}_{2}(\mathbf{k}^{\prime})+a^{s}_{2}(\mathbf{k})a^{s^{\prime}}_{1}(\mathbf{k}^{\prime})+a^{s}_{3}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle$
(49) $\displaystyle=-T_{0}\partial_{T_{3}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle+T_{0}^{\phantom{a}2}\partial_{T_{1}}\partial_{T_{2}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle-\frac{T_{0}^{\phantom{a}3}}{6}\partial_{T_{1}}^{3}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle$
$\displaystyle-T_{0}\partial_{T_{1}}\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})+a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})b^{s}_{2}(\mathbf{k})\rangle.$

From the first closure (33), the contributions proportional to $T_{0}^{2}$ and
$T_{0}^{3}$ vanish. The right-hand side of the equation involves terms of the
form:

$\partial_{T_{1}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})a^{s^{\prime\prime}}_{0}(\mathbf{k}^{\prime\prime})a^{s^{\prime\prime\prime}}_{0}(\mathbf{k}^{\prime\prime\prime})\rangle,$
(50)

which vanishes according to eq. (34). Therefore, in order to maintain a
bounded left-hand side, we must impose the asymptotic condition:

$\partial_{T_{3}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle=0,$ (51)

hence, the wave action does not evolve with respect to $T_{3}$ either. In
other words, this timescale is still too short to allow the development of a
turbulent cascade. As we will see in the next section, wave turbulence
develops on a timescale $T_{4}$.

###  III.6 Fourth asymptotic closure at time $\mathbf{T_{4}}$

At order $\epsilon^{4}$, we have the following expression:

$\displaystyle\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{4}(\mathbf{k}^{\prime})+a^{s}_{1}(\mathbf{k})a^{s^{\prime}}_{3}(\mathbf{k}^{\prime})+a^{s}_{2}(\mathbf{k})a^{s^{\prime}}_{2}(\mathbf{k}^{\prime})+a^{s}_{3}(\mathbf{k})a^{s^{\prime}}_{1}(\mathbf{k}^{\prime})+a^{s}_{4}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle$
(52) $\displaystyle=-T_{0}\partial_{T_{4}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle+\frac{T_{0}^{\phantom{a}2}}{2}\left[\partial_{T_{1}}\partial_{T_{3}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle+\partial_{T_{2}}^{2}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle\right]$
$\displaystyle-\frac{T_{0}^{\phantom{a}3}}{2}\partial_{T_{1}}^{2}\partial_{T_{2}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle+\frac{T_{0}^{\phantom{a}4}}{24}\partial_{T_{1}}^{4}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle$
$\displaystyle+\langle
b^{s}_{2}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})\rangle+\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{4}(\mathbf{k}^{\prime})+a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})b^{s}_{4}(\mathbf{k})\rangle$
$\displaystyle-T_{0}\partial_{T_{2}}\langle
a_{0}^{s}(\mathbf{k})b_{2}^{s^{\prime}}(\mathbf{k}^{\prime})+a_{0}^{s^{\prime}}(\mathbf{k}^{\prime})b_{2}^{s}(\mathbf{k})\rangle+\frac{T_{0}^{\phantom{a}2}}{2}\partial_{T_{1}}^{2}\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})+a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})b^{s}_{2}(\mathbf{k})\rangle.$

Using eq. (33), (34), (43), and (48), we can simplify the previous relation
significantly and obtain:

$\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{4}(\mathbf{k}^{\prime})+a^{s}_{1}(\mathbf{k})a^{s^{\prime}}_{3}(\mathbf{k}^{\prime})+a^{s}_{2}(\mathbf{k})a^{s^{\prime}}_{2}(\mathbf{k}^{\prime})+a^{s}_{3}(\mathbf{k})a^{s^{\prime}}_{1}(\mathbf{k}^{\prime})+a^{s}_{4}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle\\\
=-T_{0}\partial_{T_{4}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle+\langle
b^{s}_{2}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})\rangle+\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{4}(\mathbf{k}^{\prime})+a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})b^{s}_{4}(\mathbf{k})\rangle.$
(53)

Interestingly, in the expression (52), the cancellation of the term involving
the $T_{2}$ time derivative is due to the Gaussian assumption made at $t=0$.

According to the Riemann-Lebesgue lemma, both terms $\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{4}(\mathbf{k}^{\prime})+b^{s}_{4}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle$
and $\langle
b^{s}_{2}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})\rangle$ exhibit
secular divergences. Thus, the evolution of the wave action is given by:

$\partial_{T_{4}}\langle
a^{s}_{0}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle=\mathcal{C}_{T_{0}}\left\\{\langle
b^{s}_{2}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})\rangle+\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{4}(\mathbf{k}^{\prime})+a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})b^{s}_{4}(\mathbf{k})\rangle\right\\}.$
(54)

The subtlety here is that the resonance can be twofold. In the uni-resonant
case, these terms show a linear secular drift in $T_{0}$, while in the bi-
resonant case, they exhibit a drift in $T_{0}^{2}$. But it will be shown later
that:

$\mathcal{C}_{T_{0}^{2}}\left\\{\langle
b^{s}_{2}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})\rangle+\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{4}(\mathbf{k}^{\prime})+a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})b^{s}_{4}(\mathbf{k})\rangle\right\\}=0.$
(55)

The sixth-order moment $\langle
a^{s_{1}}_{0}(\mathbf{k}_{1})a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})a^{s_{4}}_{0}(\mathbf{k}_{4})a^{s_{5}}_{0}(\mathbf{k}_{5})a^{s_{6}}_{0}(\mathbf{k}_{6})\rangle$
can be written as a sum of different cumulants of the form (see Appendix A):

  * •

$q_{0}^{s_{1}s_{2}}(\mathbf{k}_{1},\mathbf{k}_{2})q_{0}^{s_{3}s_{4}}(\mathbf{k}_{3},\mathbf{k}_{4})q_{0}^{s_{5}s_{6}}(\mathbf{k}_{5},\mathbf{k}_{6})\delta^{12}(\mathbf{k})\delta^{34}(\mathbf{k})\delta^{56}(\mathbf{k})$
and other permutations;

  * •

$q_{0}^{s_{1}s_{2}}(\mathbf{k}_{1},\mathbf{k}_{2})q_{0}^{s_{3}s_{4}s_{5}s_{6}}(\mathbf{k}_{3},\mathbf{k}_{4},\mathbf{k}_{5},\mathbf{k}_{6})\delta^{12}(\mathbf{k})\delta^{3456}(\mathbf{k})$
and other permutations;

  * •

$q_{0}^{s_{1}s_{2}s_{3}}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})q_{0}^{s_{4}s_{5}s_{6}}(\mathbf{k}_{4},\mathbf{k}_{5},\mathbf{k}_{6})\delta^{123}(\mathbf{k})\delta^{456}(\mathbf{k})$
and other permutations;

  * •

$q_{0}^{s_{1}s_{2}s_{3}s_{4}s_{5}s_{6}}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4},\mathbf{k}_{5},\mathbf{k}_{6})\delta^{123456}(\mathbf{k})$.

Therefore, we need to examine all the different contributions (41 terms) to
determine which ones could lead to a secular divergence. We can immediately
disregard the contribution of
$q_{0}^{s_{1}s_{2}s_{3}s_{4}s_{5}s_{6}}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4},\mathbf{k}_{5},\mathbf{k}_{6})$
$\delta^{123456}(\mathbf{k})$ and terms of the form
$q_{0}^{s_{1}s_{2}s_{3}}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})\delta^{123}(\mathbf{k})$
$q_{0}^{s_{4}s_{5}s_{6}}(\mathbf{k}_{4},\mathbf{k}_{5},\mathbf{k}_{6})\delta^{456}(\mathbf{k})$.
This is because the different Dirac deltas do not impose any constraints on
the integration over the resonant manifold.

#### Study of the term $\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{4}(\mathbf{k}^{\prime})+b^{s}_{4}(\mathbf{k})a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})\rangle$

Using the theory of generalized functions, different types of oscillating
integrals can be evaluated, which generalizes the Riemann-Lebesgue lemma. In
particular, we can find the following long time limits of $E(X,Y,T)$ (see
Appendix C):

$E(X,Y,T)\xrightarrow[T\to+\infty]{}\begin{cases}\left[\pi\delta(X)+i\mathcal{P}\left(\frac{1}{X}\right)\right]\left[\pi\delta(Y)+i\mathcal{P}\left(\frac{1}{Y}\right)\right],&\text{if}~{}X,Y\neq
0,~{}Y\text{and}~{}X\neq Y;\\\
\left[\pi\delta(X)+i\mathcal{P}\left(\frac{1}{X}\right)\right]\left[T-i\frac{\partial}{\partial
X}\right],&\text{if}~{}X\neq 0~{}\text{and}~{}Y=0;\\\
\left[\pi\delta(Y)+i\mathcal{P}\left(\frac{1}{Y}\right)\right]\left[T-i\frac{\partial}{\partial
Y}\right],&\text{if}~{}X=0~{}\text{and}~{}Y\neq 0;\\\
\left[\pi\delta(X)+i\mathcal{P}\left(\frac{1}{X}\right)\right]i\frac{\partial}{\partial
X},&\text{if}~{}X=Y\neq 0;\\\ \frac{T^{2}}{2},&\text{if}~{}X=Y=0.\end{cases}$
(56)

We find that $\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{4}(\mathbf{k}^{\prime})\rangle$ exhibits
a secular divergence if the integration is constrained to one of the following
manifolds:

$\left\\{\begin{aligned}
s^{\prime}\omega^{\prime}&=s_{1}\omega_{1}+s_{2}\omega_{2}+s_{3}\omega_{3}\\\
s^{\prime}\omega^{\prime}&\neq
s_{2}\omega_{2}+s_{3}\omega_{3}+s_{4}\omega_{4}+s_{5}\omega_{5}+s_{6}\omega_{6};\end{aligned}\right.$
(57a) or, $\left\\{\begin{aligned} s^{\prime}\omega^{\prime}&\neq
s_{1}\omega_{1}+s_{2}\omega_{2}+s_{3}\omega_{3}\\\
s^{\prime}\omega^{\prime}&=s_{2}\omega_{2}+s_{3}\omega_{3}+s_{4}\omega_{4}+s_{5}\omega_{5}+s_{6}\omega_{6}.\end{aligned}\right.$
(57b)

There are two types of terms that satisfy the conditions (57a). The first type
is proportional to
$q_{0}^{s_{2}s_{3}}(\mathbf{k}_{2},\mathbf{k}_{3})\delta^{23}(\mathbf{k})$,
and the second type is proportional to
$q_{0}^{ss_{2}}(\mathbf{k},\mathbf{k}_{2})\delta^{02}(\mathbf{k})$ (or
$q_{0}^{ss_{3}}(\mathbf{k},\mathbf{k}3)\delta^{03}(\mathbf{k})$). It is
important to note that such a decomposition must satisfy $\Omega^{1}_{456}\neq
0$. The contribution from the first type of term is given by:

$3\delta^{00^{\prime}}(\mathbf{k})F(\mathbf{k}^{\prime},s^{\prime})\int_{\mathbb{R}^{3n}}\sum_{\begin{subarray}{c}s_{4},s_{5},s_{6}\\\
\Omega^{0^{\prime}}_{456}\neq
0\end{subarray}}\mathbf{L}^{-s^{\prime}s_{4}s_{5}s_{6}}_{-\mathbf{k}^{\prime}\mathbf{k}_{4}\mathbf{k}_{5}\mathbf{k}_{6}}\langle
a^{s}_{0}(\mathbf{k})a^{s_{4}}_{0}(\mathbf{k}_{4})a^{s_{5}}_{0}(\mathbf{k}_{5})a^{s_{6}}_{0}(\mathbf{k}_{6})\rangle\\\\[-20.0pt]
\times\delta^{0^{\prime}}_{456}(\mathbf{k})\Delta(\Omega^{0^{\prime}}_{456},T_{0})~{}\prod_{i=4}^{6}\mathrm{d}^{n}\mathbf{k}_{i}.$
(58)

On the other hand, the contribution from the second type of term
(symmetrically) is given by:

$6\langle
a_{0}^{s}(\mathbf{k})a_{0}^{s^{\prime}}(\mathbf{k}^{\prime})\rangle\int_{\mathbb{R}^{4n}}\sum_{\begin{subarray}{c}s_{1},s_{4},s_{5},s_{6}\\\
\Omega^{1}_{456}\neq
0\end{subarray}}\mathbf{L}^{-s^{\prime}s^{\prime}s_{1}-s_{1}}_{-\mathbf{k}^{\prime}\mathbf{k}^{\prime}\mathbf{k}_{1}-\mathbf{k}_{1}}\mathbf{L}^{-s_{1}s_{4}s_{5}s_{6}}_{-\mathbf{k}_{1}\mathbf{k}_{4}\mathbf{k}_{5}\mathbf{k}_{6}}\langle
a^{s_{1}}_{0}(\mathbf{k}_{1})a^{s_{4}}_{0}(\mathbf{k}_{4})a^{s_{5}}_{0}(\mathbf{k}_{5})a^{s_{6}}_{0}(\mathbf{k}_{6})\rangle\\\\[-20.0pt]
\times\delta^{1}_{456}(\mathbf{k})\Delta(\Omega^{1}_{456},T_{0})~{}\mathrm{d}^{n}\mathbf{k}_{1}\prod_{i=4}^{6}\mathrm{d}^{n}\mathbf{k}_{i}.$
(59)

A decomposition that satisfies the conditions (57b) can be expressed as:

$\langle a^{s}_{0}(\mathbf{k})a^{s_{i}}_{0}(\mathbf{k}_{i})\rangle\langle
a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{j}}_{0}(\mathbf{k}_{j})\rangle\langle
a^{s_{3}}_{0}(\mathbf{k}_{3})a^{s_{k}}_{0}(\mathbf{k}_{k})\rangle,$ (60)

where $i,j$, and $k$ are distinct elements chosen from ${4,5,6}$. It is
important to note that due to symmetry, these terms will yield the same
result, as there are six of them (corresponding to the permutations of a
three-element set). The total contribution from this term is:

$18q_{0}^{ss^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})\delta^{00^{\prime}}(\mathbf{k})\delta^{ss^{\prime}}\int_{\mathbb{R}^{3n}}\sum_{\begin{subarray}{c}s_{1},s_{2},s_{3}\\\
\Omega^{0^{\prime}}_{123}\neq
0\end{subarray}}\mathbf{L}^{-s^{\prime}s_{1}s_{2}s_{3}}_{-\mathbf{k}^{\prime}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\mathbf{L}^{-s_{1}s^{\prime}-s_{2}-s_{3}}_{-\mathbf{k}_{1}\mathbf{k}^{\prime}-\mathbf{k}_{2}-\mathbf{k}_{3}}q_{0}^{s_{2}-s_{2}}(\mathbf{k}_{2},-\mathbf{k}_{2})q_{0}^{s_{3}-s_{3}}(\mathbf{k}_{3},-\mathbf{k}_{3})\\\\[-20.0pt]
\times\delta^{0^{\prime}}_{123}(\mathbf{k})\Delta(\Omega^{0^{\prime}}_{123})~{}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}.$
(61)

In the case of a bi-resonance, $\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{4}(\mathbf{k}^{\prime})\rangle$ exhibits
a quadratic divergence in $T_{0}^{2}$, which occurs when the following
conditions are satisfied:

$\left\\{\begin{aligned}
s^{\prime}\omega^{\prime}&=s_{1}\omega_{1}+s_{2}\omega_{2}+s_{3}\omega_{3}\\\
s_{1}\omega_{1}&=s_{4}\omega_{4}+s_{5}\omega_{5}+s_{6}\omega_{6}\end{aligned}\right..$
(62)

This can only be achieved through a decomposition of the form:

$\langle a^{s}_{0}(\mathbf{k})a^{s_{i}}_{0}(\mathbf{k}_{i})\rangle\langle
a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})\rangle\langle
a^{s_{j}}_{0}(\mathbf{k}_{j})a^{s_{k}}_{0}(\mathbf{k}_{k})\rangle,$

where $i,j$, and $k$ are distinct elements chosen from ${4,5,6}$ in pairs. The
decomposition yields only three terms, which, due to symmetries, give the same
result. Furthermore, the quadratic contribution of $\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{4}(\mathbf{k}^{\prime})\rangle$ is given
by:

$\frac{9}{2}\langle
a_{0}^{s}(\mathbf{k})a_{0}^{s^{\prime}}(\mathbf{k}^{\prime})\rangle
F(\mathbf{k}^{\prime},s^{\prime})^{2}.$ (63)

#### Study of the term $\langle
b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})b^{s}_{2}(\mathbf{k})\rangle$

We will employ Riemann-Lebesgue’s lemma and the following corollary (see
Appendix C):

$\Delta(X,T)\Delta(-X,T)\xrightarrow[T\to+\infty]{}2\pi
T\delta(X)+2\mathcal{P}\left(\frac{1}{X}\right)\frac{\partial}{\partial X}.$
(64)

The term $\langle
b^{s}_{2}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})\rangle$ can
exhibit a secular drift proportional to $T_{0}$ if any of the following
conditions are satisfied:

$\left\\{\begin{aligned}
s\omega&=s_{1}\omega_{1}+s_{2}\omega_{2}+s_{3}\omega_{3}\\\
s^{\prime}\omega^{\prime}&\neq
s_{4}\omega_{4}+s_{5}\omega_{5}+s_{6}\omega_{6}\end{aligned}\right.;$ (65a)
or, $\left\\{\begin{aligned} s\omega&\neq
s_{1}\omega_{1}+s_{2}\omega_{2}+s_{3}\omega_{3}\\\
s^{\prime}\omega^{\prime}&=s_{4}\omega_{4}+s_{5}\omega_{5}+s_{6}\omega_{6}\end{aligned}\right.;$
(65b) or, $s\omega-
s_{1}\omega_{1}-s_{2}\omega_{2}-s_{3}\omega_{3}=-s^{\prime}\omega^{\prime}+s_{4}\omega_{4}+s_{5}\omega_{5}+s_{6}\omega_{6}.$
(65c)

The situation (65a) arises for decomposition of the form:

$\langle
a^{s_{i}}_{0}(\mathbf{k}_{i})a^{s_{j}}_{0}(\mathbf{k}_{j})\rangle\langle
a^{s_{k}}_{0}(\mathbf{k}_{k})a^{s_{4}}_{0}(\mathbf{k}_{4})a^{s_{5}}_{0}(\mathbf{k}_{5})a^{s_{6}}_{0}(\mathbf{k}_{6})\rangle,$
(66)

where $i$, $j$, and $k$ are elements of $\left\\{1,2,3\right\\}$ pairwise
distinct. These three decomposition yield the same result. Therefore, their
contribution is:

$3\delta^{00^{\prime}}(\mathbf{k})F(\mathbf{k},s)\int_{\mathbb{R}^{3n}}\sum_{\begin{subarray}{c}s_{4},s_{5},s_{6}\\\
\Omega^{0^{\prime}}_{456}\neq
0\end{subarray}}\mathbf{L}^{-s^{\prime}s_{4}s_{5}s_{6}}_{-\mathbf{k}^{\prime}\mathbf{k}_{4}\mathbf{k}_{5}\mathbf{k}_{6}}\langle
a^{s}_{0}(\mathbf{k})a^{s_{4}}_{0}(\mathbf{k}_{4})a^{s_{5}}_{0}(\mathbf{k}_{5})a^{s_{6}}_{0}(\mathbf{k}_{6})\rangle\\\\[-20.0pt]
\times\delta^{0^{\prime}}_{456}(\mathbf{k})\Delta(\Omega^{0^{\prime}}_{456},T_{0})~{}\prod_{i=4}^{6}\mathrm{d}^{n}\mathbf{k}_{i}.$
(67)

In the same way, we can determine the contribution corresponding to condition
(65b). This one is equal to:

$3\delta^{00^{\prime}}(\mathbf{k})F(\mathbf{k}^{\prime},s^{\prime})\int_{\mathbb{R}^{3n}}\sum_{\begin{subarray}{c}s_{1},s_{2},s_{3}\\\
\Omega^{0}_{123}\neq
0\end{subarray}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\langle
a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})a^{s_{1}}_{0}(\mathbf{k}_{1})a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})\rangle\\\\[-20.0pt]
\times\delta^{0}_{123}(\mathbf{k})\Delta(\Omega^{0}_{123},T_{0})~{}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}.$
(68)

The situation described by eq. (65c) requires a more careful treatment. It can
only occur in a decomposition of the form:

$\langle
a^{s_{i_{1}}}_{0}(\mathbf{k}_{i_{1}})a^{s_{i_{2}}}_{0}(\mathbf{k}_{i_{2}})\rangle\langle
a^{s_{j_{1}}}_{0}(\mathbf{k}_{j_{1}})a^{s_{j_{2}}}_{0}(\mathbf{k}_{j_{2}})\rangle\langle
a^{s_{k_{1}}}_{0}(\mathbf{k}_{k_{1}})a^{s_{k_{2}}}_{0}(\mathbf{k}_{k_{2}})\rangle,$
(69)

where $i_{1}$, $j_{1}$, and $k_{1}$ (respectively $i_{2}$, $j_{2}$, and
$k_{2}$) are elements of $\left\\{1,2,3\right\\}$ (respectively
$\left\\{4,5,6\right\\}$) pairwise distinct. Out of the 15 terms, only 6
satisfy this condition, and due to symmetry, they yield the same result.
Consequently, the total contribution is:

$12\pi\delta^{00^{\prime}}(\mathbf{k})\delta^{ss^{\prime}}\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\mathbf{L}^{-s^{\prime}-s_{1}-s_{2}-s_{3}}_{-\mathbf{k}^{\prime}-\mathbf{k}_{1}-\mathbf{k}_{2}-\mathbf{k}_{3}}q_{0}^{s_{1}-s_{1}}(\mathbf{k}_{1},-\mathbf{k}_{1})q_{0}^{s_{2}-s_{2}}(\mathbf{k}_{2},-\mathbf{k}_{2})\\\\[-20.0pt]
\times
q_{0}^{s_{3}-s_{3}}(\mathbf{k}_{3},-\mathbf{k}_{3})\delta(\Omega^{0}_{123})\delta^{0}_{123}(\mathbf{k})~{}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}.$
(70)

Finally, the term $\langle
b^{s}_{2}(\mathbf{k})b^{s^{\prime}}2(\mathbf{k}^{\prime})\rangle$ can exhibit
a quadratic secular drift in the case of the bi-resonance, where:

$\left\\{\begin{aligned}
s\omega&=s_{1}\omega_{1}+s_{2}\omega_{2}+s_{3}\omega_{3}\\\
s^{\prime}\omega^{\prime}&=s_{4}\omega_{4}+s_{5}\omega_{5}+s_{6}\omega_{6}\end{aligned}\right..$
(71)

This situation arises only for decomposition of the form:

$\langle
a^{s_{i_{1}}}_{0}(\mathbf{k}_{i_{1}})a^{s_{i_{2}}}_{0}(\mathbf{k}_{i_{2}})\rangle\langle
a^{s_{j_{1}}}_{0}(\mathbf{k}_{j_{1}})a^{s_{k_{1}}}_{0}(\mathbf{k}_{k_{1}})\rangle\langle
a^{s_{j_{2}}}_{0}(\mathbf{k}_{j_{2}})a^{s_{k_{2}}}_{0}(\mathbf{k}_{k_{2}})\rangle,$
(72)

where $i_{1}$, $j_{1}$, and $k_{1}$ (respectively, $i_{2}$, $j_{2}$, and
$k_{2}$) are elements of $\left\\{1,2,3\right\\}$ (respectively,
$\left\\{4,5,6\right\\}$) pairwise distinct. There are nine such
decompositions that yield the same result due to symmetry. Hence, the total
contribution of this term is given by:

$9\langle a_{0}^{s}(\mathbf{k})a_{0}^{s}(\mathbf{k}^{\prime})\rangle
F(\mathbf{k},s)F(\mathbf{k}^{\prime},s^{\prime}).$ (73)

###  III.7 Kinetic equation at time $\mathbf{T_{4}}$

All the previous results can be summarized and simplified using the symmetry
properties of the interaction term. In particular, we have:

$\mathbf{L}^{s^{\prime}s^{\prime}s_{2}-s_{2}}_{\mathbf{k}^{\prime}\mathbf{k}^{\prime}\mathbf{k}_{2}-\mathbf{k}_{2}}=-\mathbf{L}^{sss_{2}-s_{2}}_{\mathbf{k}\mathbf{k}\mathbf{k}_{2}-\mathbf{k}_{2}}.$
(74)

Hence, we finally find the following expression for the secular contributions:

$\displaystyle\mathcal{C}_{T_{0}}\left\\{\langle
b^{s}_{2}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})\rangle+\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{4}(\mathbf{k}^{\prime})+a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})b^{s}_{4}(\mathbf{k})\rangle\right\\}$
(75)
$\displaystyle=12\pi~{}\delta^{00^{\prime}}(\mathbf{k})\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\mathbf{L}^{-s^{\prime}-s_{1}-s_{2}-s_{3}}_{-\mathbf{k}^{\prime}-\mathbf{k}_{1}-\mathbf{k}_{2}-\mathbf{k}_{3}}q_{0}^{s_{1}-s_{1}}(\mathbf{k}_{1},-\mathbf{k}_{1})$
$\displaystyle\times
q_{0}^{s_{2}-s_{2}}(\mathbf{k}_{2},-\mathbf{k}_{2})q_{0}^{s_{3}-s_{3}}(\mathbf{k}_{3},-\mathbf{k}_{3})\delta(\Omega^{0}_{123})\delta^{0}_{123}(\mathbf{k})~{}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}$
$\displaystyle+18q_{0}^{ss^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})\delta^{00^{\prime}}(\mathbf{k})\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\mathbf{L}^{-s_{1}s-s_{2}-s_{3}}_{-\mathbf{k}_{1}\mathbf{k}-\mathbf{k}_{2}-\mathbf{k}_{3}}q_{0}^{s_{2}-s_{2}}(\mathbf{k}_{2},-\mathbf{k}_{2})$
$\displaystyle\times
q_{0}^{s_{3}-s_{3}}(\mathbf{k}_{3},-\mathbf{k}_{3})\delta^{0}_{123}(\mathbf{k})\Delta(\Omega^{0}_{123},T_{0})~{}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}$
$\displaystyle+18~{}q_{0}^{ss^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})\delta^{00^{\prime}}(\mathbf{k})\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-s^{\prime}s_{1}s_{2}s_{3}}_{-\mathbf{k}^{\prime}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\mathbf{L}^{-s_{1}s^{\prime}-s_{2}-s_{3}}_{-\mathbf{k}_{1}\mathbf{k}^{\prime}-\mathbf{k}_{2}-\mathbf{k}_{3}}q_{0}^{s_{2}-s_{2}}(\mathbf{k}_{2},-\mathbf{k}_{2})$
$\displaystyle\times
q_{0}^{s_{3}-s_{3}}(\mathbf{k}_{3},-\mathbf{k}_{3})\delta^{0^{\prime}}_{123}(\mathbf{k})\Delta(\Omega^{0^{\prime}}_{123},T_{0})~{}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}.$

Furthermore, we find that all the contributions in quadratic secular drifts in
$T_{0}$ vanish:

$\mathcal{C}_{T_{0}^{\phantom{a}2}}\left\\{\langle
b^{s}_{2}(\mathbf{k})b^{s^{\prime}}_{2}(\mathbf{k}^{\prime})\rangle+\langle
a^{s}_{0}(\mathbf{k})b^{s^{\prime}}_{4}(\mathbf{k}^{\prime})+a^{s^{\prime}}_{0}(\mathbf{k}^{\prime})b^{s}_{4}(\mathbf{k})\rangle\right\\}=0.$
(76)

So, the cumulants are well-ordered and the development is consistent.

Before taking the limit $T_{0}\to+\infty$, we introduce the wave action as:

$n^{s}(\mathbf{k})=q_{0}^{s-s}(\mathbf{k},-\mathbf{k}).$ (77)

Next, we integrate eq. (75) with respect to $\mathbf{k}^{\prime}$ and take the
limit $T_{0}\to+\infty$. Since the integrand decays sufficiently quickly
(thanks to the cumulant), we can interchange the limit and the integral. This
yields the following expression:

$\displaystyle\partial_{T_{4}}$ $\displaystyle n^{s}(\mathbf{k})$ (78)
$\displaystyle=12\pi\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\mathbf{L}^{s-s_{1}-s_{2}-s_{3}}_{\mathbf{k}-\mathbf{k}_{1}-\mathbf{k}_{2}-\mathbf{k}_{3}}n^{s_{1}}(\mathbf{k}_{1})n^{s_{2}}(\mathbf{k}_{2})n^{s_{3}}(\mathbf{k}_{3})\delta(\Omega^{0}_{123})\delta^{0}_{123}(\mathbf{k})~{}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}$
$\displaystyle+18\pi
n^{s}(\mathbf{k})\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{ss_{1}s_{2}s_{3}}_{\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\mathbf{L}^{-s_{1}-s-s_{2}-s_{3}}_{-\mathbf{k}_{1}-\mathbf{k}-\mathbf{k}_{2}-\mathbf{k}_{3}}n^{s_{2}}(\mathbf{k}_{2})n^{s_{3}}(\mathbf{k}_{3})\delta^{0123}(\mathbf{k})\delta(\Omega^{0123})~{}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}$
$\displaystyle+18\pi
n^{s}(\mathbf{k})\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\mathbf{L}^{-s_{1}s-s_{2}-s_{3}}_{-\mathbf{k}_{1}\mathbf{k}-\mathbf{k}_{2}-\mathbf{k}_{3}}n^{s_{2}}(\mathbf{k}_{2})n^{s_{3}}(\mathbf{k}_{3})\delta^{0}_{123}(\mathbf{k})\delta(\Omega^{0}_{123})~{}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}$
$\displaystyle=12\pi\int_{\mathbb{R}^{3n}}\sum_{s_{1},s_{2},s_{3}}\mathbf{L}^{s-s_{1}-s_{2}-s_{3}}_{\mathbf{k}-\mathbf{k}_{1}-\mathbf{k}_{2}-\mathbf{k}_{3}}\left(\frac{\mathbf{L}^{-ss_{1}s_{2}s_{3}}_{-\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}}{n^{s}(\mathbf{k})}+\frac{\mathbf{L}^{-s_{1}ss_{2}s_{3}}_{-\mathbf{k}_{1}\mathbf{k}\mathbf{k}_{2}\mathbf{k}_{3}}}{n^{s_{1}}(\mathbf{k}_{1})}+\frac{\mathbf{L}^{-s_{2}ss_{1}s_{3}}_{-\mathbf{k}_{2}\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{3}}}{n^{s_{2}}(\mathbf{k}_{2})}+\frac{\mathbf{L}^{-s_{3}ss_{1}s_{2}}_{-\mathbf{k}_{3}\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}}}{n^{s_{3}}(\mathbf{k}_{3})}\right)$
$\displaystyle\hskip 160.0pt\times
n^{s}(\mathbf{k})n^{s_{1}}(\mathbf{k}_{1})n^{s_{2}}(\mathbf{k}_{2})n^{s_{3}}(\mathbf{k}_{3})\delta(\Omega^{0}_{123})\delta^{0}_{123}(\mathbf{k})~{}\prod_{i=1}^{3}\mathrm{d}^{n}\mathbf{k}_{i}.$

The last equality has been computed by applying the changes of variables
$\mathbf{k}^{\prime}_{i}=-\mathbf{k}_{i}$ and $s^{\prime}_{i}=-s_{i}$ in the
second integral and by taking advantage of symmetries on the resonant
manifold. This last equation represents the kinetic equation for four-wave
interactions.

##  IV The case of gravitational waves

We will apply the previous formalism to gravitational waves when the Hadad-
Zakharov metric is used. Therefore, we will introduce the specific expression
of the interaction term
$\mathbf{L}^{ss_{1}s_{2}s_{3}}_{\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}$
derived earlier, into expression (78). This will allow us to investigate the
various instances of resonance in greater detail. Using the definition of
normal variables in eq. (8) it is clear that both wave action and energy do
not depend on $s$. Thus, we will simply write them as $n(\mathbf{k})$ and
$e(\mathbf{k})=\omega_{\mathbf{k}}n(\mathbf{k})$, respectively.

###  IV.1 Resonance conditions

Equation (78) provides us with the interaction conditions for a quartet of
wave vectors
$\left(\mathbf{k},\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3}\right)$. These
conditions are given by:

$\begin{cases}\mathbf{k}=\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}\\\
sk=s_{1}k_{1}+s_{2}k_{2}+s_{3}k_{3}\end{cases}.$ (79)

We can identify four possible cases:

  *   * •

All signs are the same:

$s=s_{1}=s_{2}=s_{3}.$ (80a)

  * •

Only one sign is different from the others, for instance:

$s=-s_{1}=s_{2}=s_{3}.$ (80b)

  * •

One sign is equal to $s$, and the other two signs are different, for instance:

$s=s_{1}=-s_{2}=-s_{3}.$ (80c)

  * •

All free signs are different from $s$:

$s=-s_{1}=-s_{2}=-s_{3}.$ (80d)

####  $4\leftrightarrow 0$ interactions

The $4\leftrightarrow 0$ interactions are described by the case of eq. (80d).
They are quickly ruled out. Indeed, in this situation the sum of the norms of
the wave vectors is zero, so all the wave vectors are zero. In particular
$\mathbf{k}$, which, given that
$\mathbf{L}^{ss_{1}s_{2}s_{3}}_{\mathbf{0}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}=0$,
implies that the contribution of this interaction vanishes.

####  $3\leftrightarrow 1$ interactions

The $3\leftrightarrow 1$ interactions correspond to the cases described by eq.
(80a) and (80c). In these scenarios, the wave vectors must be aligned and two
situations are possible as shown in fig. (1). The wave vectors can then be
written as follow:

$\mathbf{k}=k\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}\quad\text{and}\quad\mathbf{k}_{i}=s_{i}k_{i}\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}.$ (81)

This decomposition is introduced into the expression of
$\mathbf{L}^{ss_{1}s_{2}s_{3}}_{\mathbf{k}\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}$
to show, after a direct calculation that it vanishes. By leveraging the
different symmetries of the interaction coefficient, this result is also
proved to be true when one sign is equal to $s$ and the other two are
differents.

Figure 1: Diagram of the interaction mode in the case of eq. (80a) (on top)
and in the case of eq. (80c) (on bottom).

####  $2\leftrightarrow 2$ interactions

The $2\leftrightarrow 2$ interactions are described by the case of the
equation (80b). It allows much more geometries in the interaction quartet
because the wave vectors are not necessarily aligned, thus the wave vectors
are not simply proportional to each other. However, a change of variable in
the equation (78): $\mathbf{k}_{1}\to-\mathbf{k}_{1}$, change the resonance
condition as following:

$\begin{cases}\mathbf{k}+\mathbf{k}_{1}=\mathbf{k}_{2}+\mathbf{k}_{3}\\\
k+k_{1}=k_{2}+k_{3}\end{cases}.$ (82)

Therefore, for a given pair of wave vector
$\left(\mathbf{k},\mathbf{k}_{1}\right)$, a solution of (82) can be see as a
point on an ellipse whose focus are the extremities of $\mathbf{k}$ and
$\mathbf{k}_{1}$ (see Figure 2).

With this representation, we note that in the limit of nonlocal interactions
originating from a small wave vector (say $k\to 0$), the ellipse necessarily
tends to be flat with the foci that tend to be close to the ellipse. In this
limit, the three other wave vectors tend to be aligned, a situation in which
the interaction coefficient vanishes.

Figure 2: Schematic representation of the manifold upon $2\leftrightarrow 2$
interactions are possible. It is simply defined by the position of a dot $M$
(in red) located on an ellipse whose foci are $\mathrm{F}_{1}$ and
$\mathrm{F}_{2}$.

###  IV.2 Kinetic equation of gravitational waves

Equation (78) can be further simplified by using the following additional
symmetry, valid on the resonant manifold (82) and verified numerically:

$\mathbf{L}^{ss-
s-s}_{\mathbf{k}_{1}\mathbf{k}-\mathbf{k}_{2}-\mathbf{k}3}=\mathbf{L}^{ss-
s-s}_{\mathbf{k}\mathbf{k}_{1}-\mathbf{k}_{2}-\mathbf{k}_{3}}.$ (83)

Finally, we obtain a classical form for the kinetic equation, namely:

$\partial_{t}n(\mathbf{k})=36\pi\epsilon^{4}\int_{\mathbb{R}^{6}}\left|\mathbf{L}^{-s-sss}_{-\mathbf{k}-\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\right|^{2}\left(\frac{1}{n(\mathbf{k})}+\frac{1}{n(\mathbf{k}_{1})}-\frac{1}{n(\mathbf{k}_{2})}-\frac{1}{n(\mathbf{k}_{3})}\right)\\\
\times
n(\mathbf{k})n(\mathbf{k}_{1})n(\mathbf{k}_{2})n(\mathbf{k}_{3})\delta^{01}_{23}(\omega)\delta^{01}_{23}(\mathbf{k})\prod_{i=1}^{3}\mathrm{d}^{2}\mathbf{k}_{i}.$
(84)

We naturally use the definition $T_{4}=\epsilon^{4}t$, and the interaction
coefficient is given by:

$\displaystyle
36\pi^{2}\left|\mathbf{L}^{-s-sss}_{-\mathbf{k}-\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\right|^{2}$
(85)
$\displaystyle=\frac{\pi}{4kk_{1}k_{2}k_{3}}\Bigg{\\{}\frac{kk_{1}\left(k_{3}p_{2}q_{2}+k_{2}p_{3}q_{3}\right)-pp_{1}\left(k_{2}p_{3}q_{2}+k_{3}p_{2}q_{3}\right)-qq_{1}\left(k_{2}p_{2}q_{3}+k_{3}p_{3}q_{2}\right)}{(k_{2}+k_{3})(p_{2}+p_{3})(q_{2}+q_{3})}$
$\displaystyle+\frac{kk_{2}\left(k_{1}p_{3}q_{3}-k_{3}p_{1}q_{1}\right)+pp_{2}\left(k_{1}p_{3}q_{1}-k_{3}p_{1}q_{3}\right)+qq_{2}\left(k_{1}p_{1}q_{3}-k_{3}p_{3}q_{1}\right)}{(k_{3}-k_{1})(p_{3}-p_{1})(q_{3}-q_{1})}$
$\displaystyle+\frac{kk_{3}\left(k_{1}p_{2}q_{2}-k_{2}p_{1}q_{1}\right)+pp_{3}\left(k_{1}p_{2}q_{1}-k_{2}p_{1}q_{2}\right)+qq_{3}\left(k_{1}p_{1}q_{2}-k_{2}p_{2}q_{1}\right)}{(k_{2}-k_{1})(p_{2}-p_{1})(q_{2}-q_{1})}\Bigg{\\}}^{2}.$

The form of this interaction coefficient is not the same as that derived by
[20]. The reason is that in the Hamiltonian derivation, an additional symmetry
is introduced. However, we have verified numerically that both expressions
give the same result when this additional symmetry is also introduced in
expression (85). As often observed (as with capillary waves [45]), the
expression of the interaction coefficient obtained with the Eulerian
derivation has a more global symmetrical appearance. In the case of GW
turbulence, the difference with [20] is striking.

###  IV.3 Exact solutions

Given that this problem only allows for $2\leftrightarrow 2$ interactions, the
conservation of wave action and energy becomes apparent on equation (84).
Consequently, we can formulate two conservation equations in the spectral
domain:

$\partial_{t}n(\mathbf{k})+\bm{\nabla}_{\mathbf{k}}\Xi(\mathbf{k})=0,$ (86a)
and, $\partial_{t}e(\mathbf{k})+\bm{\nabla}_{\mathbf{k}}\Pi(\mathbf{k})=0,$
(86b)

where $\Xi(\mathbf{k})$ and $\Pi(\mathbf{k})$ are respectively the wave action
flux and the energy flux. Under the assumption of isotropic turbulence, we
integrate angularly the previous equation, we find:

$\partial_{t}N(k)+\partial_{k}\zeta(k)=0,$ (87a) and,
$\partial_{t}E(k)+\partial_{k}\varepsilon(k)=0,$ (87b)

where $N(k)=\int_{0}^{2\pi}kn(\mathbf{k})\mathrm{d}\theta=2\pi
kn(\mathbf{k})$, $E(k)=2\pi ke(\mathbf{k})$, $\zeta(k)=2\pi k\Xi(\mathbf{k})$
and $\varepsilon(k)=2\pi k\Pi(\mathbf{k})$. Furthermore, we use the isotropic
assumption in order to determine the isotropic kinetic equation by integrating
equation (84) upon the angles. We have:

$\partial_{t}N(k)=\int_{(\mathbb{R}^{+})^{3}}\mathcal{C}_{kk_{1}k_{2}k_{3}}\left(\frac{k}{N(k)}+\frac{k_{1}}{N(k_{1})}-\frac{k_{2}}{N(k_{2})}-\frac{k_{3}}{N(k_{3})}\right)\\\
\times
N(k)N(k_{1})N(k_{2})N(k_{3})\delta^{01}_{23}(\omega)\mathrm{d}k_{1}\mathrm{d}k_{2}\mathrm{d}k_{3},$
(88)

where:

$\mathcal{C}_{kk_{1}k_{2}k_{3}}=\frac{9\epsilon^{4}}{\pi}\int_{0}^{2\pi}\int_{0}^{2\pi}\int_{0}^{2\pi}\int_{0}^{2\pi}\left|\mathbf{L}^{-s-sss}_{-\mathbf{k}-\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}\right|^{2}\delta^{01}_{23}(\mathbf{k})\mathrm{d}\theta\mathrm{d}\theta_{1}\mathrm{d}\theta_{2}\mathrm{d}\theta_{3}.$
(89)

A subtler point to note here is that the integration in equation (88) is
constrained by Dirac’s delta function $\delta^{01}_{23}(\omega)$ and the
definition of the isotropic coefficient in equation (89). Therefore, it’s not
just an integration over $(\mathbb{R}^{+})^{3}$, but rather an integration
over a more complex manifold. It’s important to remember that for
gravitational waves in our unit system, $\omega=k$, so
$\delta^{01}_{23}(\omega)$ becomes $\delta^{01}_{23}(k)$.

We look for scale-invariant solutions so that, we assume: $N(k)=Ak^{x}$ where
$A>0$ and $x$ are two real constants. We also introduce the dimensionless
wavenumber $\xi_{i}=k_{i}/k$. Equation (88) can thus be written as follow:

$\partial_{t}N(k)=A^{3}k^{3x+1}\int_{(\mathbb{R}^{+})^{3}}\mathcal{C}_{1\xi_{1}\xi_{2}\xi_{3}}\left(\xi_{1}\xi_{2}\xi_{3}\right)^{x}\left(1+\xi_{1}^{\phantom{a}1-x}-\xi_{2}^{\phantom{a}1-x}-\xi_{3}^{\phantom{a}1-x}\right)\\\\[-20.0pt]
\times\delta(1+\xi_{1}-\xi_{2}-\xi_{3})\mathrm{d}\xi_{1}\mathrm{d}\xi_{2}\mathrm{d}\xi_{3}.$
(90)

Following [20], it is clear that if:

$x=0\quad\text{or}\quad x=1,$ (91)

the integral vanishes. These solutions correspond to thermodynamic equilibrium
in which both the wave action and energy fluxes are zero. Nonetheless, they
are not the only accessible solutions. We perform a Zakharov’s transform in
order to determine other solutions. We will duplicate four times the integral
in the right hand side of the equation (90) and we will perform the following
variable changes leaving one of them unchanged:

$\displaystyle\mathrm{TZa:}\left(\xi_{1},\xi_{2},\xi_{3}\right)$
$\displaystyle\longrightarrow\left(\frac{1}{\xi_{1}},\frac{\xi_{2}}{\xi_{1}},\frac{\xi_{3}}{\xi_{1}}\right);$
(92a) $\displaystyle\mathrm{TZb:}\left(\xi_{1},\xi_{2},\xi_{3}\right)$
$\displaystyle\longrightarrow\left(\frac{\xi_{3}}{\xi_{2}},\frac{1}{\xi_{2}},\frac{\xi_{1}}{\xi_{2}}\right);$
(92b) $\displaystyle\mathrm{TZc:}\left(\xi_{1},\xi_{2},\xi_{3}\right)$
$\displaystyle\longrightarrow\left(\frac{\xi_{2}}{\xi_{3}},\frac{\xi_{1}}{\xi_{3}},\frac{1}{\xi_{3}}\right).$
(92c)

We finally have:

$\partial_{t}N(k)=\frac{A^{3}}{4}k^{3x+1}I(x),$ (93)

where:

$I(x)=\int_{(\mathbb{R}^{+})^{3}}\mathcal{C}_{1\xi_{1}\xi_{2}\xi_{3}}\left(\xi_{1}\xi_{2}\xi_{3}\right)^{x}\left(1+\xi_{1}^{\phantom{a}1-x}-\xi_{2}^{\phantom{a}1-x}-\xi_{3}^{\phantom{a}1-x}\right)\\\\[-10.0pt]
\times\left(1+\xi_{1}^{\phantom{a}-3x-2}-\xi_{2}^{\phantom{a}-3x-2}-\xi_{3}^{\phantom{a}-3x-2}\right)\delta(1+\xi_{1}-\xi_{2}-\xi_{3})\mathrm{d}\xi_{1}\mathrm{d}\xi_{2}\mathrm{d}\xi_{3}.$
(94)

Therefore we find two other solutions:

$x=-2/3\quad\text{or}\quad x=-1.$ (95)

They have non-zero constant fluxes which indicates the occurrence of cascades.
To determine the direction of the cascades, we need to use equations (87a) and
(87b). For the case $x=-2/3$, the flux can be expressed as:

$\partial_{t}N(k)=\frac{A^{3}}{4}k^{3x+1}I(x)=-\partial_{k}\zeta(k),$ (96)

so that:

$\zeta(k)=-\frac{A^{3}}{4}\frac{k^{3x+2}}{3x+2}I(x).$ (97)

Therefore, in this case, we have:

$\zeta=-\frac{A^{3}}{4}\int_{(\mathbb{R}^{+})^{3}}\mathcal{C}_{1\xi_{1}\xi_{2}\xi_{3}}\left(\xi_{1}\xi_{2}\xi_{3}\right)^{-2/3}\left(1+\xi_{1}^{5/3}-\xi_{2}^{5/3}-\xi_{3}^{5/3}\right)\ln\left(\frac{\xi_{2}\xi_{3}}{\xi_{1}}\right)\\\
\times\delta\left(1+\xi_{1}-\xi_{2}-\xi_{3}\right)\mathrm{d}\xi_{1}\mathrm{d}\xi_{2}\mathrm{d}\xi_{3},$
(98)

which is negative as a numerical evaluation shows. This indicates the presence
of an inverse cascade. Similarly, for the case $x=-1$, we have:

$\partial_{t}E(k)=\frac{A^{3}}{4}k^{3x+2}I(x)=-\partial_{k}\varepsilon(k)$
(99)

so that:

$\varepsilon(k)=-\frac{A^{3}}{4}\frac{k^{3x+3}}{3x+3}I(x),$ (100)

which, at $x=-1$, yields the energy flow:

$\varepsilon=-\frac{A^{3}}{4}\int_{(\mathbb{R}^{+})^{3}}\mathcal{C}_{1\xi_{1}\xi_{2}\xi_{3}}\left(\xi_{1}\xi_{2}\xi_{3}\right)^{-1}\left(1+\xi_{1}^{2}-\xi_{2}^{2}-\xi_{3}^{2}\right)\left(-\xi_{1}\ln\xi_{1}+\xi_{2}\ln\xi_{2}+\xi_{3}\ln\xi_{3}\right)\\\
\times\delta\left(1+\xi_{1}-\xi_{2}-\xi_{3}\right)\mathrm{d}\xi_{1}\mathrm{d}\xi_{2}\mathrm{d}\xi_{3},$
(101)

which is positive as shown by a numerical evaluation, indicating a direct
cascade.

In conclusion, the equations of the system (7d) give rise to two cascades with
respective Kolmogorov-Zakharov spectra of the form:

$N(k)\propto\left(-\zeta\right)^{1/3}k^{-2/3}\quad\text{and}\quad
N(k)\propto\varepsilon^{1/3}k^{-1}.$ (102)

##  V Conclusion

Wave turbulence is a successful analytical theory with many applications [33,
34]. Our main objective was to use the multiple time scale method originally
proposed by Benney and Newell [22] to derive an integro-differential equation,
known as the kinetic equation, for quartic gravity (surface) wave
interactions. Unlike the seminal paper, our derivation is very general and
assumes only a few symmetries for the interaction coefficient. Consequently,
our result can be applied to any problem where nonlinearities are cubic and
turbulence statistically homogeneous.

A key aspect of this approach is to identify secular drifts arising from the
different decompositions of spectral cumulants. It is these drifts that enable
us to establish consistent long-term dynamics and derive the wave kinetic
equation. In practice, this requires the evaluation of oscillating integrals
via generalized functions theory. We have shown that a natural closure is
obtained because problematic terms such as sixth-order cumulants are not
secular. Therefore, in the long-time limit (which also means $\epsilon\to 0$),
these contributions become asymptotically small at main order. Note that this
conclusion does not exclude the possibility of sixth-order cumulants
contributing to higher orders in the development.

Furthermore, we have shown that quartic wave turbulence retains a memory up to
time $T_{3}$ in the sense that the evolution of the fourth-order cumulant at
zero order depends on the initial condition, so that if it is initially zero,
it reaches time $T_{3}$ without evolving. This property can even be
generalized to higher order.

By applying the multiple time scale method to gravitational waves, we have
obtained the kinetic equation that describes the temporal evolution of a set
of waves of weak amplitude. The dynamics of gravitational waves is slow,
involving a typical time scale proportional to $\tau_{GW}/\epsilon^{4}$, where
$\tau_{GW}\sim 1/\omega$ is a linear time and $\epsilon$ is a small parameter
measuring the amplitude of gravitational waves. The kinetic equation
asymptotically describes the transfer of wave action and energy through
Fourier modes. Exact stationary solutions for isotropic turbulence can be
obtained using the Zakharov transformation. These solutions are called
Kolmogorov-Zakharov spectra. For gravitational wave turbulence, they
correspond to a direct cascade of energy and an inverse cascade of wave
action. As far as we know, these are the first exact solutions of general
relativity of a statistical nature.

Our study was initially based on an Eulerian derivation of the wave amplitude
equation. This makes a difference with the first study based on a Hamiltonian
derivation [20]. The kinetic equations obtained are slightly different in
their expression, yet retain the same degree of homogeneity. Consequently, the
exact solutions are the same, as is the direction of the cascades. The
difference is understood as the consequence of an additional symmetry
introduced in the Hamiltonian derivation. Both approaches are therefore fully
compatible.

Direct numerical proof of the existence of a dual cascade in gravitational
wave turbulence has recently been obtained [43]. For the future, it is
important to continue this numerical study in order to verify the power law
indices and to measure, if possible, the acceleration of the inverse cascade
[41]. Another topic concerns intermittency. In the language of turbulence,
this means, on the one hand, checking how far wave amplitudes are from a
Gaussian distribution (a small difference is expected) and, on the other,
measuring structure functions. The latter involves measuring the difference
between a field taken at two positions separated by a distance $L$, all at a
given power, and observing the dependence of these functions on $L$. Power
laws are expected. These basic fields could be the components of the space-
time metric.

###### Acknowledgements.

We acknowledge Sergey Nazarenko for useful discussions.

##  Appendix A Expression of the six points cummulants

Assuming the means field is equal to zero, the sixth-order moment can be
decomposed on the spectral cumulants according to:

$\displaystyle\langle
a^{s_{1}}_{0}(\mathbf{k}_{1})a^{s_{2}}_{0}(\mathbf{k}_{2})a^{s_{3}}_{0}(\mathbf{k}_{3})a^{s_{4}}_{0}(\mathbf{k}_{4})a^{s_{5}}_{0}(\mathbf{k}_{5})a^{s_{6}}_{0}(\mathbf{k}_{6})\rangle$
(103)
$\displaystyle=q_{0}^{12}(\mathbf{k})q_{0}^{34}(\mathbf{k})q_{0}^{56}(\mathbf{k})+q_{0}^{12}(\mathbf{k})q_{0}^{35}(\mathbf{k})q_{0}^{46}(\mathbf{k})+q_{0}^{12}(\mathbf{k})q_{0}^{36}(\mathbf{k})q_{0}^{45}(\mathbf{k})+q_{0}^{13}(\mathbf{k})q_{0}^{24}(\mathbf{k})q_{0}^{56}(\mathbf{k})$
$\displaystyle+q_{0}^{13}(\mathbf{k})q_{0}^{25}(\mathbf{k})q_{0}^{46}(\mathbf{k})+q_{0}^{13}(\mathbf{k})q_{0}^{26}(\mathbf{k})q_{0}^{45}(\mathbf{k})+q_{0}^{14}(\mathbf{k})q_{0}^{23}(\mathbf{k})q_{0}^{56}(\mathbf{k})+q_{0}^{14}(\mathbf{k})q_{0}^{25}(\mathbf{k})q_{0}^{36}(\mathbf{k})$
$\displaystyle+q_{0}^{14}(\mathbf{k})q_{0}^{26}(\mathbf{k})q_{0}^{35}(\mathbf{k})+q_{0}^{15}(\mathbf{k})q_{0}^{23}(\mathbf{k})q_{0}^{46}(\mathbf{k})+q_{0}^{15}(\mathbf{k})q_{0}^{24}(\mathbf{k})q_{0}^{36}(\mathbf{k})+q_{0}^{15}(\mathbf{k})q_{0}^{26}(\mathbf{k})q_{0}^{34}(\mathbf{k})$
$\displaystyle+q_{0}^{16}(\mathbf{k})q_{0}^{23}(\mathbf{k})q_{0}^{45}(\mathbf{k})+q_{0}^{16}(\mathbf{k})q_{0}^{24}(\mathbf{k})q_{0}^{35}(\mathbf{k})+q_{0}^{16}(\mathbf{k})q_{0}^{25}(\mathbf{k})q_{0}^{34}(\mathbf{k})$
$\displaystyle+q_{0}^{123}(\mathbf{k})q_{0}^{456}(\mathbf{k})+q_{0}^{124}(\mathbf{k})q_{0}^{356}(\mathbf{k})+q_{0}^{125}(\mathbf{k})q_{0}^{346}(\mathbf{k})+q_{0}^{126}(\mathbf{k})q_{0}^{345}(\mathbf{k})+q_{0}^{134}(\mathbf{k})q_{0}^{256}(\mathbf{k})$
$\displaystyle+q_{0}^{135}(\mathbf{k})q_{0}^{246}(\mathbf{k})+q_{0}^{136}(\mathbf{k})q_{0}^{246}(\mathbf{k})+q_{0}^{145}(\mathbf{k})q_{0}^{236}(\mathbf{k})+q_{0}^{146}(\mathbf{k})q_{0}^{235}(\mathbf{k})+q_{0}^{156}(\mathbf{k})q_{0}^{234}(\mathbf{k})$
$\displaystyle+q_{0}^{12}(\mathbf{k})q_{0}^{3456}(\mathbf{k})+q_{0}^{13}(\mathbf{k})q_{0}^{2456}(\mathbf{k})+q_{0}^{14}(\mathbf{k})q_{0}^{2356}(\mathbf{k})+q_{0}^{15}(\mathbf{k})q_{0}^{2346}(\mathbf{k})+q_{0}^{16}(\mathbf{k})q_{0}^{2345}(\mathbf{k})$
$\displaystyle+q_{0}^{23}(\mathbf{k})q_{0}^{1456}(\mathbf{k})+q_{0}^{24}(\mathbf{k})q_{0}^{1356}(\mathbf{k})+q_{0}^{25}(\mathbf{k})q_{0}^{1346}(\mathbf{k})+q_{0}^{26}(\mathbf{k})q_{0}^{1345}(\mathbf{k})+q_{0}^{34}(\mathbf{k})q_{0}^{1256}(\mathbf{k})$
$\displaystyle+q_{0}^{35}(\mathbf{k})q_{0}^{1246}(\mathbf{k})+q_{0}^{36}(\mathbf{k})q_{0}^{1245}(\mathbf{k})+q_{0}^{45}(\mathbf{k})q_{0}^{1236}(\mathbf{k})+q_{0}^{46}(\mathbf{k})q_{0}^{1235}(\mathbf{k})+q_{0}^{56}(\mathbf{k})q_{0}^{1234}(\mathbf{k})$
$\displaystyle+q_{0}^{123456}(\mathbf{k}),$

where $q_{0}^{i_{1}i_{2}\dots
i_{r}}(\mathbf{k})=q_{0}^{s_{i_{1}}s_{i_{2}}\dots
s_{i_{r}}}(\mathbf{k}_{i_{1}},\mathbf{k}_{i_{2}},\dots,\mathbf{k}_{i_{r}})\delta\left(\sum_{j=1}^{r}\mathbf{k}_{i_{j}}\right)$.

##  Appendix B Riemann-Lebesgue’s lemma for distributions

We define for $(x,t)\in\mathbb{R}\times\mathbb{R}^{+}$,
$\Delta(x,t)=\int_{0}^{t}e^{ixt^{\prime}}~{}\mathrm{d}t^{\prime}=\frac{e^{ixt}-1}{ix}$.
In terms of generalized function, $\Delta(\cdot,t)$ has the following
asymptotic behaviour:

$\Delta(\cdot,t)\xrightarrow[t\to+\infty]{}\begin{cases}\pi\delta(\cdot)+i\mathcal{P}\left(\frac{1}{\cdot}\right),&\text{if}~{}x\neq
0;\\\ t,&\text{if}~{}x=0.\\\ \end{cases}$ (104)

Proof \- Let $\varphi:x\mapsto\varphi(x)$ be a test function of
$C_{c}^{\infty}(\mathbb{R})$, this mean that $\varphi$ is infinitely
continuously differentiable with a compact support. Thus, we define
$M\in\mathbb{R}$ such that:

$\forall|x|\geq M,~{}\varphi(x)=0.$ (105)

We also have:

$\displaystyle\langle\Delta(\cdot,t),\varphi\rangle$
$\displaystyle=\int_{\mathbb{R}}\frac{e^{ixt}-1}{ix}\varphi(x)\mathrm{d}x$
(106)
$\displaystyle=\int_{\mathbb{R}}\left[\frac{\sin(xt)}{x}\varphi(x)+i\frac{1-\cos(xt)}{x}\varphi(x)\right]\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}}\frac{\sin(u)}{u}\varphi\left(\frac{u}{t}\right)\mathrm{d}u+i\int_{0}^{+\infty}\frac{1-\cos(xt)}{x}\left[\varphi(x)-\varphi(-x)\right]\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}}\frac{\sin(u)}{u}\varphi\left(\frac{u}{t}\right)\mathrm{d}u+i\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(x)}{x}\mathrm{d}x\right\\}-i\int_{0}^{+\infty}\frac{\cos(xt)}{x}\left[\varphi(x)-\varphi(-x)\right]~{}\mathrm{d}x.$

Concerning the real part, it is clear that:

$\forall
x\in\mathbb{R},~{}\frac{\sin(u)}{u}\varphi\left(\frac{u}{t}\right)\xrightarrow[t\to+\infty]{}\frac{\sin(u)}{u}\varphi\left(0\right),$
(107)

and:

$\forall
x\in\mathbb{R},~{}\left|\frac{\sin(u)}{u}\varphi\left(\frac{u}{t}\right)\right|\leq\begin{cases}0,&\text{
if }|x|\geq M;\\\ \max\limits_{\mathbb{R}}|\varphi|,&\text{ else};\end{cases}$
(108)

which is integrable. Thus, according to the dominated convergence theorem, the
real part converges to:

$\varphi(0)\int_{\mathbb{R}}\frac{\sin(u)}{u}~{}\mathrm{d}u=\pi\varphi(0).$
(109)

Concerning the imaginary part, we may apply the Riemann-Lebesgue’s theorem for
functions since $v\mapsto\frac{\varphi(v)-\varphi(-v)}{v}$ is integrable upon
$\mathbb{R}^{+}$. Hence,

$\int_{0}^{+\infty}\frac{\cos(vt)}{v}\left[\varphi(v)-\varphi(-v)\right]~{}\mathrm{d}v\xrightarrow[t\to+\infty]{}0.$
(110)

By gathering all these elements, we have:

$\langle\Delta(\cdot,t),\varphi\rangle\xrightarrow[t\to+\infty]{}\pi\varphi(0)+i\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(v)}{v}\mathrm{d}v\right\\}.$
(111)

Furthermore, it is straightforward to see that if $x=0$, $\Delta(x,t)=t$.
$\square$

##  Appendix C Consequences of the lemma: other asymptotic properties

We now define for $(x,y,t)\in\mathbb{R}^{2}\times\mathbb{R}^{+}$,
$E\left(x,y,t\right)=\int_{0}^{t}\Delta(x-y,t^{\prime})e^{iyt^{\prime}}~{}\mathrm{d}t^{\prime}$.
In terms of distribution, $E(\cdot,\cdot,t)$ has the following asymptotic
behaviour:

$E(x,y,t)\xrightarrow[t\to+\infty]{}\begin{cases}\left[\pi\delta(x)+i\mathcal{P}\left(\frac{1}{x}\right)\right]\left[\pi\delta(y)+i\mathcal{P}\left(\frac{1}{y}\right)\right],&\text{if}~{}x\neq
0,~{}y\neq 0~{}\text{and}~{}x\neq y;\\\
\left[\pi\delta(x)+i\mathcal{P}\left(\frac{1}{x}\right)\right]\left[t-i\frac{\partial}{\partial
x}\right],&\text{if}~{}x\neq 0~{}\text{and}~{}y=0;\\\
i\left[\pi\delta(x)+i\mathcal{P}\left(\frac{1}{x}\right)\right]\frac{\partial}{\partial
x},&\text{if}~{}x=y\neq 0;\\\ \frac{t^{2}}{2},&\text{if}~{}x=y=0.\end{cases}$
(112)

Proof \- Let $\varphi:(x,y)\mapsto\varphi(x,y)$ be a test function of
$C_{c}^{\infty}(\mathbb{R}^{2},\mathbb{R})$.
$\bullet$ We first assume that both $x$ and $y$ are not constrained to vanish
and are not constrained to be equal. We have:

$\displaystyle\langle E(\cdot,\cdot,t),\varphi\rangle$
$\displaystyle=\int_{\mathbb{R}^{2}}\frac{\Delta(x,t)-\Delta(y,t)}{i(x-y)}\varphi(x,y)\mathrm{d}x\mathrm{d}y$
(113)
$\displaystyle=\int_{\mathbb{R}}\Delta(x,t)\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(x,y)}{i(x-y)}\mathrm{d}y\right\\}\mathrm{d}x-\int_{\mathbb{R}}\Delta(y,t)\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(x,y)}{i(x-y)}\mathrm{d}x\right\\}\mathrm{d}y.$

Thus, when taking the limit $t\to+\infty$, we obtain:

$\displaystyle\langle
E(\cdot,\cdot,t),\varphi\rangle\xrightarrow[t\to+\infty]{}$
$\displaystyle~{}i\pi\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(0,y)}{y}\mathrm{d}y\right\\}+\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{1}{x}\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(x,y)}{(x-y)}\mathrm{d}y\right\\}\mathrm{d}x\right\\}$
(114)
$\displaystyle+i\pi\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(x,0)}{x}\mathrm{d}x\right\\}+\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{1}{y}\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(x,y)}{(x-y)}\mathrm{d}x\right\\}\mathrm{d}y\right\\}.$

Hence, when using the Poincaré–Bertrand lemma, we have:

$\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{1}{x}\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(x,y)}{(x-y)}\mathrm{d}y\right\\}\mathrm{d}x\right\\}+\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{1}{y}\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(x,y)}{(x-y)}\mathrm{d}x\right\\}\mathrm{d}y\right\\}\\\
=-\pi^{2}\varphi(0,0)+\mathcal{P}\left\\{\int_{\mathbb{R}^{2}}\frac{\varphi(x,y)}{xy}\mathrm{d}x\mathrm{d}y\right\\}.$
(115)

This result allows us to write:

$\langle
E(\cdot,\cdot,t),\varphi\rangle\xrightarrow[t\to+\infty]{}~{}-\pi^{2}\varphi(0,0)+\mathcal{P}\left\\{\int_{\mathbb{R}^{2}}\frac{\varphi(x,y)}{xy}\mathrm{d}x\mathrm{d}y\right\\}\\\
+i\pi\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(0,y)}{y}\mathrm{d}y\right\\}+i\pi\mathcal{P}\left\\{\int_{\mathbb{R}}\frac{\varphi(x,0)}{x}\mathrm{d}x\right\\}.$
(116)

Finally, if $x\neq 0$, $y\neq 0$ and $x\neq y$, we have:

$E(\cdot,\cdot,t)\xrightarrow[t\to+\infty]{}\left[\pi\delta(x)+i\mathcal{P}\left(\frac{1}{x}\right)\right]\left[\pi\delta(y)+i\mathcal{P}\left(\frac{1}{y}\right)\right].$
(117)

$\bullet$ Now, we assume that only $y$ is constrained to be equal to zero. In
this case, we have:

$\displaystyle\langle E(\cdot,0,t),\varphi\rangle$
$\displaystyle=\int_{\mathbb{R}}\frac{e^{ixt}-1-ixt}{(ix)^{2}}\varphi(x)\mathrm{d}x$
(118)
$\displaystyle=\int_{\mathbb{R}}-\frac{1}{x}\left[e^{ixt}\frac{\mathrm{d}\varphi}{\mathrm{d}x}(x)+it\varphi(x)e^{ixt}\right]+\frac{1}{x}\frac{\mathrm{d}\varphi}{\mathrm{d}x}(x)-\frac{t}{ix}\varphi(x)\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}}\Delta(x,t)\left(t-i\frac{\mathrm{d}}{\mathrm{d}x}\right)\varphi(x)\mathrm{d}x.$

Thus, when taking the limit $t\to+\infty$ we find the asymptotic behaviour:

$E(\cdot,0,t)\xrightarrow[t\to+\infty]{}\left[\pi\delta(x)+i\mathcal{P}\left(\frac{1}{x}\right)\right]\left(t-i\frac{\mathrm{d}}{\mathrm{d}x}\right).$
(119)

$\bullet$ When $y$ is equal to $x$, we may notice that:

$E(x,x,t)=t\Delta(x,t)-E(x,0,t).$ (120)

Thus, the asymptotic behaviour is given by:

$E(x,x,t)\xrightarrow[t\to+\infty]{}i\left[\pi\delta(x)+i\mathcal{P}\left(\frac{1}{x}\right)\right]\frac{\mathrm{d}}{\mathrm{d}x}.$
(121)

$\bullet$ Last but not least, it is clear that if $x=y=0$, we have
$E(\cdot,\cdot,t)=t^{2}/2$. $\square$

In terms of generalized function, $\Delta(x,t)\Delta(-x,t)$ has the following
asymptotic behaviour:

$\Delta(x,t)\Delta(-x,t)\xrightarrow[t\to+\infty]{}2\pi
t\delta(x)+2\mathcal{P}\left(\frac{1}{x}\right)\frac{\mathrm{d}}{\mathrm{d}x}.$
(122)

Proof \- We may notice that:

$\Delta(x,t)\Delta(-x,t)=E(x,0,t)+E(-x,0,t).$ (123)

Yet, we have:

$E(-x,0,t)\xrightarrow{t\to+\infty}\left[\pi\delta(x)-i\mathcal{P}\left(\frac{1}{x}\right)\right]\left(t+i\frac{\mathrm{d}}{\mathrm{d}x}\right).$
(124)

Hence, the asymptotic behaviour is given by:

$\Delta(x,t)\Delta(-x,t)\xrightarrow{t\to+\infty}2\pi
t\delta(x)+2\mathcal{P}\left(\frac{1}{x}\right)\frac{\mathrm{d}}{\mathrm{d}x}.\quad\square$
(125)

## References

  * Newell [1968] A. C. Newell, Reviews of Geophysics 6, 1 (1968).
  * LIGO and Virgo Scientific Collaborations [2016] LIGO and Virgo Scientific Collaborations, Physical Review Letters 116, 061102 (2016).
  * Einstein [1918] A. Einstein, Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften , 154 (1918).
  * Maggiore [2007] M. Maggiore,  _Gravitational Waves: Volume 1: Theory and Experiments_ (OUP Oxford, Oxford, 2007).
  * Krauss [1992] L. M. Krauss, Physics Letters B 284, 229 (1992).
  * Kosowsky  _et al._ [1992] A. Kosowsky, M. S. Turner, and R. Watkins, Physical Review Letters 69, 2026 (1992).
  * Kamionkowski  _et al._ [1994] M. Kamionkowski, A. Kosowsky, and M. S. Turner, Physical Review D 49, 2837 (1994).
  * Fenu  _et al._ [2009] E. Fenu, D. G. Figueroa, R. Durrer, and J. García-Bellido, Journal of Cosmology and Astroparticle Physics 2009 (10), 005.
  * Rubakov  _et al._ [1982] V. A. Rubakov, M. V. Sazhin, and A. V. Veryaskin, Physics Letters B 115, 189 (1982).
  * Guzzetti  _et al._ [2016] M. C. Guzzetti, N. Bartolo, M. Liguori, and S. Matarrese, La Rivista del Nuovo Cimento 39, 399 (2016), arxiv:1605.01615 [astro-ph] .
  * Collaboration [2016] P. Collaboration, Astronomy & Astrophysics 594, A20 (2016), arxiv:1502.02114 [astro-ph] .
  * Bertone  _et al._ [2023] G. Bertone, O. L. Buchmueller, and P. S. Cole, npj Microgravity 9, 1 (2023).
  * Collaboration [2023] N. Collaboration, The Astrophysical Journal 951, L8 (2023).
  * Zakharov [1967] V. E. Zakharov, Soviet Journal of Experimental and Theoretical Physics 24, 455 (1967).
  * Galtier [2003] S. Galtier, Physical Review E 68, 015301(R) (2003).
  * David and Galtier [2022] V. David and S. Galtier, Journal of Plasma Physics 88, 905880509 (2022).
  * Düring  _et al._ [2006] G. Düring, C. Josserand, and S. Rica, Physical Review Letters 97, 025503 (2006).
  * Nazarenko and Onorato [2006] S. Nazarenko and M. Onorato, Physica D: Nonlinear Phenomena 219, 1 (2006), arxiv:nlin/0507051 .
  * Dyachenko  _et al._ [1992] S. Dyachenko, A. C. Newell, A. Pushkarev, and V. E. Zakharov, Physica D: Nonlinear Phenomena 57, 96 (1992).
  * Galtier and Nazarenko [2017] S. Galtier and S. V. Nazarenko, Physical Review Letters 119, 221101 (2017).
  * Kraichnan [1959] R. H. Kraichnan, Journal of Fluid Mechanics 5, 497 (1959).
  * Benney  _et al._ [1966] D. J. Benney, P. G. Saffman, and G. K. Batchelor, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 289, 301 (1966).
  * Poincaré [1892] H. Poincaré,  _Les méthodes nouvelles de la mécanique céleste (Vol. 2)_ (Villars-Gauthier Editeur, 1892).
  * Benney and Newell [1967] D. J. Benney and A. C. Newell, The Physics of Fluids 10, S281 (1967).
  * Sagdeev and Galeev [1969] R. Z. Sagdeev and A. Galeev,  _Nonlinear Plasma Theory_ (W. A. Benjamin, 1969).
  * Vedenov [1967] A. A. Vedenov, in  _Reviews of Plasma Physics: Volume 3_, edited by M. A. Leontovich (Springer US, Boston, MA, 1967) pp. 229–276.
  * Zakharov [1965] V. E. Zakharov, Journal of Applied Mechanics and Technical Physics 6, 22 (1965).
  * Zakharov and Filonenko [1967a] V. E. Zakharov and N. N. Filonenko, Journal of Applied Mechanics and Technical Physics 8, 37 (1967a).
  * Zakharov and Filonenko [1967b] V. Zakharov and N. N. Filonenko (1967).
  * Kaner and Yakovenko [1970] E. A. Kaner and V. M. Yakovenko, Zh. Eksp. Teor. Fiz. 58: 587-96(Feb 1970). (1970).
  * Bizoń and Rostworowski [2011] P. Bizoń and A. Rostworowski, Physical Review Letters 107, 031102 (2011), arxiv:1104.3702 [gr-qc, physics:hep-th] .
  * de Oliveira  _et al._ [2013] H. P. de Oliveira, L. A. Pando Zayas, and E. L. Rodrigues, Physical Review Letters 111, 051101 (2013).
  * Nazarenko [2011] S. Nazarenko,  _Wave Turbulence_, Lecture Notes in Physics, Vol. 825 (Springer, Berlin, Heidelberg, 2011).
  * Galtier [2022] S. Galtier,  _Physics of Wave Turbulence_ (Cambridge University Press, Cambridge, 2022).
  * Dias  _et al._ [2012a] O. J. C. Dias, G. T. Horowitz, D. Marolf, and J. E. Santos, Classical and Quantum Gravity 29, 235019 (2012a), arxiv:1208.5772 [gr-qc, physics:hep-th] .
  * Dias  _et al._ [2012b] O. J. C. Dias, G. T. Horowitz, and J. E. Santos, Classical and Quantum Gravity 29, 194002 (2012b), arxiv:1109.1825 [gr-qc, physics:hep-th] .
  * Craps  _et al._ [2014] B. Craps, O. Evnin, and J. Vanhoof, Journal of High Energy Physics 2014, 48 (2014), arxiv:1407.6273 [gr-qc, physics:hep-th, physics:nlin] .
  * Craps  _et al._ [2015] B. Craps, O. Evnin, and J. Vanhoof, Journal of High Energy Physics 2015, 108 (2015), arxiv:1412.3249 [gr-qc, physics:hep-th, physics:nlin] .
  * Hadad and Zakharov [2014] Y. Hadad and V. Zakharov, Journal of Geometry and Physics 80, 37 (2014), arxiv:1401.7251 [gr-qc] .
  * Galtier  _et al._ [2019] S. Galtier, S. V. Nazarenko, É. Buchlin, and S. Thalabard, Physica D: Nonlinear Phenomena 390, 84 (2019).
  * Galtier  _et al._ [2020] S. Galtier, J. Laurie, and S. V. Nazarenko, Universe 6, 98 (2020).
  * Collaboration [2020] P. Collaboration, Astronomy & Astrophysics 641, A10 (2020).
  * Galtier and Nazarenko [2021] S. Galtier and S. V. Nazarenko, Physical Review Letters 127, 131101 (2021).
  * Deng and Hani [2022] Y. Deng and Z. Hani, Propagation of chaos and the higher order statistics in the wave kinetic theory (2022), arxiv:2110.04565 [math-ph] .
  * Galtier [2021] S. Galtier, Geophysical and Astrophysical Fluid Dynamics 115, 234 (2021).