Basile Gallet

Sergey Nazarenko

Bérengère Dubrulle

Wave-turbulence description of interacting particles: Klein-Gordon model with a Mexican-hat potential

Keywords: number(s): 05, 45, -a, 11, 10, Lm

In field theory, particles are waves or excitations that propagate on the fundamental state. In experiments or cosmological models, one typically wants to compute the out-of-equilibrium evolution of a given initial distribution of such waves. Wave turbulence deals with out-of-equilibrium ensembles of weakly nonlinear waves, and is therefore well suited to address this problem. As an example, we consider the complex Klein-Gordon equation with a Mexican-hat potential. This simple equation displays two kinds of excitations around the fundamental state: massive particles and massless Goldstone bosons. The former are waves with a nonzero frequency for vanishing wave number, whereas the latter obey an acoustic dispersion relation. Using waveturbulence theory, we derive wave kinetic equations that govern the coupled evolution of the spectra of massive and massless waves. We first consider the thermodynamic solutions to these equations and study the wave condensation transition, which is the classical equivalent of Bose-Einstein condensation. We then focus on nonlocal interactions in wave-number space: we study the decay of an ensemble of massive particles into massless ones. Under rather general conditions, these massless particles accumulate at low wave number. We study the dynamics of waves coexisting with such a strong condensate, and we compute rigorously a nonlocal Kolmogorov-Zakharov solution, where particles are transferred nonlocally to the condensate, while energy cascades towards large wave numbers through local interactions. This nonlocal cascading state constitutes the intermediate asymptotics between the initial distribution of waves and the thermodynamic state reached in the long-time limit.

I. INTRODUCTION

Some nonlinear partial differential equations (PDEs) are simple enough and yet they describe deep and nontrivial universal processes in a broad range of physical applications. Well-known examples of such PDEs are the Korteweg-de Vries and the Gross-Pitaevskii (GP) equations. It would not be an exaggeration to put into such a class a universal model that allows to understand the concept of spontaneous breaking of global symmetry. Consider, for instance, the following complex nonlinear Klein-Gordon equation

ψ tt -ψ + (-1 + |ψ| 2 )ψ = 0 , (1)
where ψ(x,t) ∈ C. The linear part of this equation resembles the standard Klein-Gordon equation, where the square mass coefficient m 2 is replaced by -1. Such a negative square mass (m 2 = -1) was originally believed to lead to superluminal particles. However, it was soon realized that such a field cannot sustain a localized particlelike excitation because the base state ψ = 0 is unstable. Equation (1) conserves the total energy

E = |ψ t | 2 + |∇ψ| 2 -|ψ| 2 + 1 2 |ψ| 4 dx , ( 2
)
where the last two terms correspond to the famous Mexicanhat potential U (ψ) = -|ψ| 2 + |ψ| 4 /2; we therefore refer to Eq. ( 1) as the Klein-Gordon Mexican-hat (KGMH) model. Equation ( 1) is invariant to a uniform shift in phase of ψ, i.e., ψ → e iα ψ, with α = const. The invariant corresponding to this U(1) symmetry is the "charge"

Q = i 2 (ψ∂ t ψ * -ψ * ∂ t ψ) dx , (3)
where ψ * denotes the complex conjugate of ψ. In this paper, we focus mostly on the situation Q = 0, Q = 0 being briefly described in Sec. III C. In high energy physics, the KGMH model is a nonlinear σ model [1,2] with a Mexican-hat potential. It is a simple relativistic model that describes the interaction between two fields: the first one corresponds to a massive particle, while the second one corresponds to a massless Goldstone boson. As such, the KGMH model provides a simplified description of a system of sigma mesons interacting with pions: the mass of the pions is neglected in this model, as it is much less than the one of the sigma meson.

In this paper, we study the KGMH system within the framework of wave turbulence, with direct numerical simulations (DNSs) to check the predictions. Wave turbulence [START_REF] Falkovich | Kolmogorov Spectra of Turbulence[END_REF][START_REF] Nazarenko | Wave Turbulence[END_REF] deals with ensembles of weakly interacting waves. It proceeds through the derivation of a kinetic equation that governs the time evolution of the spectral density of the wave field. Among its solutions are thermal equilibrium and generalized Rayleigh-Jeans distributions. However, the main use of the kinetic equation is to describe out-of-equilibrium situations. Indeed, for scale-invariant systems, one can compute outof-equilibrium power-law solutions to the kinetic equation: these are the celebrated Kolmogorov-Zakharov spectra. More generally, the kinetic equation governs the evolution of an ensemble of waves (or particles), starting from an arbitrary weakly nonlinear initial condition.

Such out-of-equilibrium initial conditions are encountered in particle physics and cosmology: they correspond for instance to heavy-ion collisions produced experimentally, or to the early universe. There has therefore been a significant recent interest in applying the wave-turbulence approach in these domains [5,6].

From the point of view of wave-turbulence theory, the KGMH model is a new object of study, which is interesting and important for several reasons: spontaneous symmetry breaking leads to two kinds of waves, or particles, with distinct dispersion relations. This contrasts with most wave-turbulence systems, where one usually deals with a single kind of waves, as, e.g., in the GP equation. In addition to transfers of energy from one scale to another, the KGMH system displays "reactions": decay of massive particles into massless ones, and fusion of massless particles into massive ones. This brings features to the wave-turbulence dynamics, with strongly nonlocal interactions between different kinds of particles, that prohibit standard Kolmogorov-Zakharov cascades.

Nevertheless, the KGMH model shares some similarities with the simpler GP equation: both have two positive invariants. Indeed, in addition to energy, the KGMH equation has an adiabatic invariant which is similar to the number of particles. Here, "adiabatic" means that this invariant is conserved exactly in the asymptotic limit of small nonlinearity only. This adiabatic invariance is sufficient to trigger wave condensation [7], a phenomenon that we study both theoretically and numerically. In the particle physics context, this corresponds to Bose-Einstein condensation of pions [2,8,9].

The decay of the massive mode (σ meson) into massless Goldstone bosons (pions) is the leading thread of this paper. In Sec. II, we propose a simple mechanical analog of the KGMH model (1), which provides physical insight into the underlying dynamics. In Sec. III, we study the linear and weakly nonlinear dynamics around the minimum of the Mexican-hat potential (the uniform fundamental state). We compute the dispersion relations for the massive and massless modes, together with a reduced Hamiltonian which takes into account the dominant three-wave interaction between these modes. This dominant interaction involves two massless modes and one massive mode. We analyze it in Sec. IV by deriving a simple set of ordinary differential equations (ODEs) for a single triad of interacting waves, which highlights the decay instability of a monochromatic massive wave into two massless monochromatic waves.

In Sec. V, we develop the wave-turbulence description of the KGMH model. We compute kinetic equations governing the evolution of the spectra of massive and massless modes. We identify two invariants of these equations: the energy, and an adiabatic invariant somewhat similar to a number of particles. The kinetic equations admit Rayleigh-Jeans equilibrium solutions, but no local Kolmogorov-Zakharov outof-equilibrium solutions. We focus on such thermodynamic equilibria in Sec. VI. For low enough energy per particle, the thermodynamic state displays wave condensation, the classical counterpart to the Bose-Einstein condensation (BEC) of pions. These findings are confirmed by DNSs, which show unambiguously that BEC is observed over a very long time, even in systems where the number of particles is only an adiabatic invariant.

The out-of-equilibrium dynamics involve strongly nonlocal processes which we consider in Secs. VII and VIII. The first section addresses the decay instability for isotropic ensembles of massive waves: on the one hand, the kinetic equations allow to predict the distribution of the massless waves that appear during the linear stage of the instability. On the other hand, the thermodynamics of Sec. VI predict the the reaction yield of the decay instability: starting from massive waves (σ mesons) at a single wave number, we predict the fraction of these waves that has decayed into massless waves (pions) in the long-time limit.

Section VIII finally addresses the behavior of ensembles of waves coexisting with a strong condensate of massless particles. We derive rigorously a reduced kinetic equation that describes a nonlocal transfer of particles to the condensate, together with a local Kolmogorov-Zakharov energy cascade towards high wave numbers. We see this cascading state as an intermediate asymptotics that describes the dynamical process of condensation: it makes the link between the initial distribution of particles and the final thermodynamic condensed state.

II. AN ANALOGOUS MECHANICAL SYSTEM

One can gain some insight into the physics described by Eq. ( 1) by considering a system from classical mechanics that follows the same equation. An example of such system is the square lattice of pendulums sketched in Fig. 1. The pendulums can move in two angular directions ζ and ξ . The corresponding kinetic energy of a pendulum is J ( ζ 2 + ξ 2 )/2, with J the moment of inertia of a pendulum. Each pendulum consists of a magnetized bead attached to a rigid rod, the magnetic dipole of the bead being aligned with the rod. Neighboring pendulums interact through two mechanisms: first, two dipoles repulse each other. Second, we assume that the fixation points of the pendulums are connected by torsional springs. Dipole-dipole repulsion efficiently propagates longitudinal perturbations through the lattice, while torsional FIG. 1. (Color online) A lattice of magnetic pendulums: small magnets placed under the vertical position of each pendulum create the Mexican-hat potential. The interaction between neighboring pendulums results from dipole-dipole interaction and torsional springs between their fixation points. In the continuous limit, this system follows the KGMH equation (1).

springs propagate transverse oscillations. The stiffness C of the torsional springs is chosen to match the strength of the dipole-dipole interaction, so that the interaction energy between the pendulums reads as

C 2 (i,j ) neighbors (ζ i -ζ j ) 2 + (ξ i -ξ j ) 2 , ( 4
)
where ζ n and ξ n are the angular degrees of freedom of the nth pendulum, and we consider dipole-dipole interactions between nearest neighbors only. We further assume that the relative displacements of the pendulums are small and consider the lowest order description of the dipole-dipole interaction. Finally, to create the Mexican-hat potential, we place a magnet under the rest position of each bead. Gravity drives a given pendulum towards the vertical, but the repulsion between the bead and the magnet displaces the equilibrium position away from the vertical, at a distance r = r 0 from the vertical axis. The combined effects of gravity and of the magnet result in a Mexican-hat-shaped potential V (r). If the length of the rods is much longer than r 0 , we can use the approximate relation r 2 2 (ζ 2 + ξ 2 ). In the neighborhood of r 0 , we approximate the potential by

V (r) -V 2 r 2 + V 4 r 4 = -V 2 2 (ζ 2 + ξ 2 ) + V 4 4 (ζ 2 + ξ 2 ) 2 .
The Lagrangian of the system follows from the subtraction of the Mexican-hat potential energy and neighboring interaction term to the kinetic term

L = n J 2 ζn 2 + ξn 2 -V 4 4 ζ 2 n + ξ 2 n 2 + V 2 2 ζ 2 n + ξ 2 n - C 2 (i,j ) neighbors (ζ i -ζ j ) 2 + (ξ i -ξ j ) 2 . ( 5
)
We consider the continuous limit and introduce the complex variable ψ(x,y,t) = ζ (x,y,t) + iξ (x,y,t). After a rescaling of time, space, and ψ, this Lagrangian takes the following dimensionless form:

L = |ψ t | 2 -|∇ψ| 2 + |ψ| 2 - 1 2 |ψ| 4 d 2 x , (6)
which corresponds to Eq. ( 1) in two spatial dimensions. In this analogous system, the charge Q corresponds to the sum of the vertical angular momenta of each pendulum around their fixation points. A nonzero value of Q means that the initial condition (and the subsequent evolution) favors one direction of azimuthal rotation of the pendulums.

III. WEAKLY NONLINEAR DYNAMICS AROUND THE MINIMUM OF THE POTENTIAL

A key feature of the KGMH model is that it displays spontaneous U(1) symmetry breaking: the fundamental state, defined as the state of lowest energy, is not the usual ψ = 0. Instead, it is a uniform state ψ = ψ 0 = const, with |ψ 0 | = 1: ψ is uniform, time independent, and lies in the minimum of the Mexican-hat potential. We wish to study the dynamics of the system around this fundamental state.

Up to a phase shift in the definition of ψ, we can choose ψ 0 = 1 and consider perturbations of this state Substitution into Eq. ( 1) yields

ψ(x,t) = 1 + φ(x,t) .
φ tt -φ + φ + φ * + 2φφ * + φ 2 + |φ| 2 φ = 0 . (7)
We focus on the weakly nonlinear dynamics arising from the quadratic terms of this equation, and therefore neglect the cubic nonlinearity. In variables λ = Re(φ) and χ = Im(φ) we have

λ tt -λ + 2λ + 3λ 2 + χ 2 = 0 , ( 8a
)
χ tt -χ + 2λχ = 0 . ( 8b
)
The linear parts of these two equations provide two dispersion relations. For a plane wave proportional to e i(ωt-k•x) , they give, respectively,

ω a k = 2 + k 2 , ( 9a
)
ω b k = k, ( 9b
)
where k = |k|. We refer to these two dispersion relations as branch a and branch b. They are plotted in Fig. 2. Branch a has a nonzero frequency for vanishing wave number, ω a = 0 for k = 0. In terms of particle physics, this mode has a nonzero rest energy and thus corresponds to massive particles. By contrast, the frequency of b waves vanishes for k = 0: branch b corresponds to a massless Goldstone boson. This is a basic model for coexisting sigma meson (branch a) and pion (branch b) fields, in which the pion's mass is neglected.

In the lattice of pendulums sketched in Fig. 1, the massive mode corresponds to oscillations of the pendulums along their radial direction. In this mode, each pendulum feels two restoring forces, one from the potential well V (r) and one from the coupling to its neighbors. By contrast, mode b corresponds to oscillations of the pendulums along their azimuthal directions. The restoring force originates only from the coupling with the neighboring pendulums, and if all the pendulums move in phase along their azimuthal direction (b mode with infinite wavelength), then there is no restoring force: mode b is massless.

A. Hamiltonian formulation

The system (8a) and (8b) follows from the Hamiltonian

H = 1 2 p 2 + 1 2 q 2 + 1 2 (∇λ) 2 + λ 2 + 1 2 (∇χ ) 2 + λ 3 + λχ 2 dx , ( 10
)
where p and q are the momenta conjugate to the variables λ and χ , respectively. Hamilton equations are

λ t = δ H δp , p t = - δ H δλ , ( 11a
)
χ t = δ H δq , q t = - δ H δχ . ( 11b
)
We consider these equations in a d-dimensional periodic cube of side L and introduce Fourier series

⎛ ⎜ ⎜ ⎜ ⎝ λ(x,t) χ (x,t) p(x,t) q(x,t) ⎞ ⎟ ⎟ ⎟ ⎠ = k ⎛ ⎜ ⎜ ⎜ ⎝ λk (t) χk (t) pk (t) qk (t) ⎞ ⎟ ⎟ ⎟ ⎠ e ik•x , ( 12
)
where the sums are over k

∈ 2π L Z d . Because these four fields are real valued, ( λ-k , χ-k , p-k , q-k ) = ( λ * k , χ * k , p * k , q * k ). In Fourier space, Hamilton equations therefore become λk = δH δ p * k , ṗk = - δH δ λ * k , ( 13a
) χk = δH δ q * k , qk = - δH δ χ * k , ( 13b
)
with

H = H 2 + H int , (14)
where

H 2 = 1 2 k | pk | 2 + ω a k 2 | λk | 2 + | qk | 2 + ω b k 2 | χk | 2
(15) and

H int = k 1 ,k 2 ,k 3 λk 1 λk 2 λk 3 + λk 1 χk 2 χk 3 δ(k 1 + k 2 + k 3 ), ( 16
)
where δ is a Kroenecker symbol, i.e., δ(k = 0) = 1 and δ(k = 0) = 0.

B. Normal variables

Let us introduce the normal variables a k and b k that diagonalize H 2 :

λk = a k + a * -k 2ω a k , pk = ω a k (a k -a * -k ) i √ 2 , (17a) χk = b k + b * -k 2ω b k , qk = ω b k (b k -b * -k ) i √ 2 . (17b)
We refer to a k and b k as the respective amplitudes of a and b waves (or modes, or particles).

In these variables the equations of motion are

i ȧk = δH δa * k , ( 18a
) i ḃk = δH δb * k , ( 18b
)
where H = H 2 + H int , with

H 2 = k ω a k |a k | 2 + ω b k |b k | 2 (19)
and

H int = 1 2 √ 2 k 1 ,k 2 ,k 3 1 ω a 1 ω a 2 ω a 3 [a 1 a 2 a 3 δ(k 1 + k 2 + k 3 ) + 3a 1 a 2 a * 3 δ(k 1 + k 2 -k 3 )] + 1 2 √ 2 k 1 ,k 2 ,k 3 1 ω a 1 ω b 2 ω b 3 [a 1 b 2 b 3 δ(k 1 + k 2 + k 3 ) + 2a 1 b 2 b * 3 δ(k 1 + k 2 -k 3 ) + a * 1 b 2 b 3 δ(k 1 -k 2 -k 3 )] + c.c., (20)
where c.c. stands for "complex conjugate," and we introduced the shorthand notations

a i = a k i , b i = b k i , ω a i = ω a k i , and
ω b i = ω b k i .
The terms in the Hamiltonian (20) correspond, respectively, to three-wave processes of type a + a + a → 0, a + a → a, a + b + b → 0, a + b → b, and b + b → a. In a weakly nonlinear description of the system, the only relevant nonlinearities are those for which both a wave number and a frequency resonance condition can be satisfied. The other nonlinearities produce nonresonant terms only, which do not contribute to the weakly nonlinear dynamics. We show in Appendix A that, out of the processes outlined above, the only one allowed by the frequency resonance condition is b + b → a. We can therefore discard all the other terms in H int , as they do not contribute to the weakly nonlinear dynamics, and simply write

H int = k 1 ,k 2 ,k 3 V 1 23 a * 1 b 2 b 3 δ 3 12 + c.c., (21)
where we have switched to a shorthand notation of Kroenecker deltas

δ 1 23 ≡ δ(k 1 -k 2 -k 3 )
, and introduced the interaction coefficient

V 1 23 ≡ 1 2 √ 2 1 ω a 1 ω b 2 ω b 3 . ( 22
)
The equations of motion become

ȧk = -iω a k a k -i k 1 ,k 2 V k 12 b 1 b 2 δ k 12 , ( 23a
) ḃk = -iω b k b k -2i k 1 ,k 2 V 1 2k a 1 b * 2 δ 1 2k . ( 23b
)
A similar Hamiltonian system with two types of interacting waves was studied in [10]. However, in that study the dominant process was of the type a + b → b, which leads to very different dynamics than the process b + b → a considered here.

C. Finite charge Q = 0

We have considered so far only perturbations around the minimum energy state, which is uniform ψ lying in the minimum of the Mexican-hat potential. This corresponds to a vanishing charge invariant Q = 0. We now briefly consider nonzero values of Q. Writing the field in polar coordinates ψ = r(x,t)e iθ(x,t) , the charge reads as

Q = r 2 θ dx . ( 24
)
A nonzero charge therefore indicates a preferred direction of rotation in the complex plane. Uniform states with nonzero charge can be sought in the form ψ 0 (t) = Re iω 0 t , where R and ω 0 are real constants. Substitution into the KGMH equation ( 1) leads to

R 2 = 1 + ω 2 0 , (25)
i.e., there is a continuous family of uniform solutions to the KGMH equation, corresponding to different values of the charge invariant

Q = R 2 ω 0 L d .
Once again, the lattice of magnetic pendulums provides a simple interpretation for these states. From Eq. (24) the charge is the sum of the vertical angular momenta of the pendulums about their fixation points. The uniform solution with ω 0 = 0 corresponds to all the pendulums parallel and at rest in the minimum of the Mexican-hat potential, whereas uniform solutions with ω 0 = 0 correspond to all the pendulums parallel and rotating at the same angular frequency ω 0 around the vertical. For the latter solutions R 2 = 1 + ω 2 0 > 1: the field is not exactly in the minimum of the Mexican-hat potential, but instead it is "bobsleighing" on its outer edge. In the system of pendulums, this "bobsleighing" results from the centrifugal force acting on each pendulum.

To study perturbations around the uniform state with nonzero Q, we write

ψ = ψ 0 (t)[1 + φ(x,t)] = √ 1 + ω 2 0 e iω 0 t [1 + φ(x,t)],
and consider infinitesimal perturbations φ 1. Substitution into the KGMH equation gives at linear order in φ

φ tt + 2iω 0 φ t -φ + 1 + ω 2 0 (φ + φ * ) = 0, ( 26
)
and after taking the real and imaginary parts, with φ = λ + iχ,

λ tt -λ + 2 1 + ω 2 0 λ -2ω 0 χ t = 0, ( 27a
)
χ tt -χ + 2ω 0 λ t = 0. ( 27b
)
As compared to the situation Q = 0, the fields λ and χ are now coupled at linear order by terms proportional to ω 0 . Inserting a plane wave structure for these two fields and asking for nontrivial solutions yields the following dispersion relation:

ω 4 -2ω 2 k 2 + 1 + 3ω 2 0 + k 2 k 2 + 2 + 2ω 2 0 = 0, (28) which has two branches of solutions ω 2 ± = k 2 + 1 + 3ω 2 0 ± 1 + 3ω 2 0 2 + 4k 2 ω 2 0 . ( 29
)
In the limit ω 0 → 0, we see that ω + branch corresponds to a modes, whereas ω -corresponds to b modes. For ω 0 = 0, the mode ω + remains massive, with lim k→0 ω + = 0, and the ω -branch remains massless, with lim k→0 ω -= 0. Note that the speed of the low-k b particles is now reduced, ∂ω -/∂k = √ (1 + ω 2 0 )/(1 + 3ω 2 0 ), and the a particle's mass is

√ 1 + 3ω 2 0
times greater than in the ω 0 = 0 case. At the linear level, the structure of the problem is nevertheless quite similar to the Q = 0 situation.

Interesting phenomena appear at the nonlinear level. As for the Q = 0 situation, one can compute normal variables a k and b k diagonalizing the quadratic part of the Hamiltonian a k (resp. b k ) corresponding to ω + (resp. ω -). These normal variables are now superpositions of the λ and χ fields, together with the corresponding momenta:

λk = ω 2 + -k 2 √ 2ω + a k + a * -k 2 - ω 2 --k 2 2ω 3 -ω 2 0 b k -b * -k 2i , (30a) μk = ω 0 √ ω + a k -a * -k i √ 2 + 1 √ ω - b k + b * -k √ 2 , ( 30b
) pk = √ ω + (ω 2 + -k 2 ) √ 2 a k -a * -k 2i + ω 2 --k 2 2ω -ω 2 0 b k + b * -k 2 , ( 30c
) qk = -ω 0 ω 3 - a k + a * -k √ 2 + √ ω - b k -b * -k i √ 2 . ( 30d
)
The weakly nonlinear limit is still governed by three-wave interactions, but the dispersion relations (29) now allow for new processes: a + a → a and a + b → b are still impossible, but, in addition to b + b → a, the processes a → a + b and a + b → a are compatible with the dispersion relation. This possibly triggers interesting new dynamics. However, because the normal variables a k and b k have rather intricate expressions when Q = 0, the computation of the kinetic equation is more involved than for Q = 0. The wave-turbulence analysis of this paper therefore focuses on Q = 0, and we leave the case Q = 0 for future work.

IV. DECAY OF THE MASSIVE MODE INTO TWO GOLDSTONE WAVES

The dominant three-wave process of the weakly nonlinear KGMH dynamics can be highlighted by studying a single triad of interacting modes. Consider the dynamical equations in the weak-nonlinearity limit, (8a) and (8b), and expand the fields as

λ = λ 0 (t,T ) + 2 λ 1 (t,T ) + O( 3 ), ( 31a
) χ = χ 0 (t,T ) + 2 χ 1 (t,T ) + O( 3 ), ( 31b
)
where the slow time scale is T = t. To order , Eqs. (8a) and (8b) reduce to

∂ tt λ 0 -λ 0 + 2λ 0 = 0, ( 32a
)
∂ tt χ 0 -χ 0 = 0. ( 32b
)
In Appendix A, we saw that triads of interacting mode are made of two b modes and one a mode. We consider a single such triad in the solution to the order O( ) equations:

λ 0 = A(T )e i( √ k 2 +2 t-k•x) + c.c., (33a)
χ 0 = B 1 (T )e i(k 1 t-k 1 •x) + B 2 (T )e i(k 2 t-k 2 •x) + c.c. (33b)
To obtain nontrivial dynamics at next order, we further assume that the resonance conditions

√ k 2 + 2 = k 1 + k 2 and k = k 1 + k 2 are
met. To order 2 , the equations become

∂ tt λ 1 -λ 1 + 2λ 1 = -2∂ tT λ 0 -3λ 2 0 -χ 2 0 , (34a) ∂ tt χ 1 -χ 1 = -2∂ tT χ 0 -2λ 0 χ 0 . ( 34b
)
The solvability condition requires the right-hand side of these equations to have no terms that are resonant with the operator on the left-hand side. We insert solutions (33a) and (33b) into Eqs. (34a) and (34b), before demanding that the right-hand side of (34a) has no term proportional to e i( √ k 2 +2 t-k•x) , and that the right-hand side of (34b) has no terms proportional to e i(k 1 t-k 1 •x) or to e i(k 2 t-k 2 •x) . These three constraints lead to the following system of ODEs for the slow evolution of the wave amplitudes:

d T A = i 2 √ k 2 + 2 B 1 B 2 , ( 35a
)
d T B 1 = i K 1 AB * 2 , ( 35b
)
d T B 2 = i K 2 AB * 1 . ( 35c
)
Consider a single b wave: (A = B 2 = 0, B 1 = 0) is a time-independent solution to the system of ODEs. We study the stability of this b wave by considering infinitesimal perturbations A 1 and B 2 1. The linearized evolution equation for these perturbations is

d T T F + |B 1 | 2 2k 2 √ k 2 + 2 F = 0, (36)
where F = A or B 2 . The b wave is stable, in the sense that the perturbations of A and B 2 oscillate but remain small. By contrast, if one starts off with a single a wave (A = 0, B 1 = B 2 = 0) and perturbs it with infinitesimal B 1 and B 2 , the subsequent evolution of these perturbations is governed by

d T T B 1;2 - |A| 2 k 1 k 2 B 1;2 = 0, ( 37
)
which admits some exponentially growing solutions: the a wave is unstable and decays spontaneously into two b waves.

In terms of particles physics, this instability would describe the disintegration of a massive particle (e.g., the σ meson) into two massless Goldstone bosons (two pions).

The growth rate γ = |A|/ √ k 1 k 2 of this instability is a first indication of nonlocal interactions in Fourier space: recalling the resonance conditions, γ is maximal when

k 1 k 2 = 1/2 [see condition (A2)], i.e., when k 1 = ( √ 2 + k 2 + k)/2 and k 2 = ( √ 2 + k 2 -k)/2 (or vice versa). Nonlocal interactions occur in at least two cases: if k 1, then k 1 k 2 √ 2 k. If k 1, then either k 1 k or k 2 k. Such nonlocal interactions
are an important feature of wave turbulence in the KGMH model.

To study the subsequent evolution of this instability, we solve Eq. ( 1) numerically. We focus on two-dimensional (2D) numerical simulations, which correspond to the dimension of the equivalent lattice of pendulums. We truncate the system at large wave numbers. Such a truncation naturally arises from the discreteness of the lattice of pendulums. For more general applications of Eq. ( 1), the truncation can be seen as a crude modelization of the ultraviolet cutoff arising from Bose statistics.

Note that the domain size is an important parameter of this problem: the dispersion relation ω a k is not self-similar, with nonrelativistic waves for k 1 and relativistic ones for k 1. To accurately describe a mixture of relativistic and nonrelativistic waves, we need to choose a domain size that is much larger than 2π . The numerical code is described in Appendix D: it uses a standard pseudospectral method, with an exact integration of the stiff linear wave operator and Adams-Bashforth time stepping to deal with the nonlinear terms. Dealiasing is performed using the one-half rule for this equation with cubic nonlinearity.

In this section, the domain size is [0,20π ] 2 with a resolution of 256 2 . After dealiasing, we keep 128 modes in each direction: k ∈ 1 10 × {-64; 64} 2 . In Fig. 3, we show snapshots of the λ and χ fields, together with their one-dimensional (1D) spectra, i.e., their spectra integrated over the direction of k. The initial condition is a small-scale progressive a wave with k = (0; 3), plus some weak background noise. The instability rapidly sets in and b waves appear. As can be seen in Eq. (37), the maximum growth rate of the instability is achieved when the product k 1 k 2 of the two b wave numbers is small. Minimizing this quantity while satisfying the resonance conditions suggests that k 1 1 and k 2 ∼ k. Indeed, in Fig. 3 we observe the emergence a combination of large-scale b waves, together with b waves around the wave number k of the initial a wave.

By contrast, in Fig. 4 we show similar snapshots for a simulation initiated with a large-scale a wave, with k = (0; 0.3). Once again, b waves rapidly appear as a result of the decay instability. However, for such a large-scale a wave, the frequency resonance condition imposes that one of the b waves has a frequency (and therefore a wave number) of order unity. But, then the resonance of the three wave vectors imposes that the second b wave also has a wave number of order unity and almost antiparallel to that of the first b wave. As a consequence, in Fig. 4 the instability sets in with b waves having wave numbers of order unity.

In the following, we extend the description of these nonlocal interactions to random ensembles of waves, using the framework of wave turbulence.

V. WAVE TURBULENCE

We now consider ensembles of weakly nonlinear waves with homogeneous statistics in an infinite domain. For such systems, wave turbulence can be used to derive wave kinetic equations that govern the evolution of the wave spectra.

A. Wave kinetic equations

Let us define the wave action spectra of the modes a and b:

n a k = L 2π d |a k | 2 , n b k = L 2π d |b k | 2 ,
where the angular brackets denote an ensemble average over many realizations.

In Appendix B, we apply standard wave-turbulence techniques to Eqs. (23a) and (23b), which leads to the following kinetic equations for the evolution of the two spectra:

ṅa k = 4π V k 12 2 n b 1 n b 2 -2n b 1 n a k δ ω abb k12 δ k 12 dk 1 dk 2 , ( 38a
)
ṅb k = 8π V 1 2k 2 n a 1 n b 2 + n a 1 n b k -n b 2 n b k δ ω abb 12k δ 1 2k dk 1 dk 2 .
(38b)

B. Conservation laws for the kinetic equations

Kinetic equations (38a) and (38b) conserve the wave energy

E = L d ω a k n a k + ω b k n b k dk, ( 39
)
as well as the following particle invariant:

N = L d 2n a k + n b k dk. ( 40
)
This second invariant is nontrivial, and it can be guessed only once the dominant three-wave process is identified. Indeed, the process a b + b is similar to a chemical reaction where a is a "molecule" made of two b "atoms," and N would be the total number of such b "atoms": one per b wave, and two per a wave. N is also the number of b particles that would be in the system if the reaction a → b + b were complete. For brevity, in the following we simply refer to N as the number of particles.

An important comment to make here is that N is only an adiabatic invariant: it is conserved by the weakly nonlinear dynamics, but there is no similar invariant for the strongly nonlinear KGMH equation. By contrast, E is the quadratic part of the global invariant E, the latter being conserved even in the strongly nonlinear regime.

C. Thermodynamic equilibrium

One can easily check that the kinetic equations admit solutions corresponding to energy equipartition in k space,

n a k = T ω a k , n b k = T ω b k , ( 41
)
for any temperature T = const, 1 and solutions corresponding to an equipartition of the particle invariant

n b k = 2n a k = const. ( 42
)
More generally, the equilibrium Rayleigh-Jeans (or thermodynamic) solution is given by

n a k = T ω a k + 2μ , n b k = T ω b k + μ , ( 43
)
where μ = const is a chemical potential. Finally, we remark that any initial condition with a vanishing spectrum n b k = 0 and an arbitrary spectrum n a k is a stationary solution to the kinetic equations (38a) and (38b). One can also show that spectra with b waves at low wave number only are also steady solutions to the kinetic equations. Such distributions with only a single type of particles are described in details in Secs. VII and VIII. and δ(ω abb k12 ) ≡ δ(ω abb k12 ) are independent of the directions of k 1 and k 2 , all that needs to be integrated over the angular variables is δ k 12 . This leads to

ṅa k = 4π D a (k) S 12k δ ω abb k12 ω a k ω b 1 ω b 2 n b 1 n b 2 -2n b 1 n a k (k 1 k 2 ) d-1 dk 1 dk 2 , ( 44a
)
ṅb k = 8π D b (k) S 12k δ ω abb 12k ω a 1 ω b k ω b 2 n a 1 n b 2 + n a 1 n b k -n b 2 n b k × (k 1 k 2 ) d-1 dk 1 dk 2 , ( 44b
)
where

S 12k = 1 2 2 k 2 k 2 1 + k 2 k 2 2 +k 2 1 k 2 2 -k 4 -k 4 1 -k 4 2 for d = 2 (45)
and

S 12k = 1 2kk 1 k 2 for d = 3. ( 46
)
The integration domains D a (k and k 2 and the additional restriction of type (A2):

D a (k) = (k 1 ,k 2 ) ∈ (R + ) 2 |k 1 k 2 1 2 ; k 1 + k 2 k; k + k 1 k 2 ; k + k 2 k 1 , ( 47a
)
D b (k) = (k 1 ,k 2 ) ∈ (R + ) 2 | k 2 1 2k ; k 1 + k 2 k; k + k 1 k 2 ; k + k 2 k 1 . ( 47b
)
They are sketched, respectively, in Figs. 5 and6.

Substituting

ω b 1 ω b 2 = k 1 k 2 and ω b k ω b 2 = kk 2 , we have ṅa k = 4π ω a k D a (k) S 12k δ ω abb k12 n b 1 n b 2 -2n b 1 n a k (k 1 k 2 ) d-2 dk 1 dk 2 , ( 48a
) ṅb k = 8π k D b (k) k 1 ω a 1 S 12k δ ω abb 12k n a 1 n b 2 + n a 1 n b k -n b 2 n b k × (k 1 k 2 ) d-2 dk 1 dk 2 . ( 48b
)
Substitution of S 12k for d = 3 leads to an even simpler system

ṅa k = 2π kω a k D a (k) δ ω abb k12 n b 1 n b 2 -2n b 1 n a k dk 1 dk 2 , ( 49a
) ṅb k = 4π k 2 D b (k) k 1 ω a 1 δ ω abb 12k n a 1 n b 2 + n a 1 n b k -n b 2 n b k dk 1 dk 2 . ( 49b

E. Absence of Kolmogorov-Zakharov spectra

The central feature of wave-turbulence theory is that one can compute power-law spectra satisfying the kinetic equations. Such Kolmogorov-Zakharov (KZ) spectra are obtained for scale-invariant systems, i.e., when the dispersion relations and interaction coefficients are homogeneous functions of their arguments. In the KGMH model, this is the case only in the relativistic limit k 1. Let us focus on this limit and exclude 2D systems: wave-turbulence kinetic equations are known to be ill posed in 2D for nondispersive waves ω ∼ k (see e.g. [11]). In 3D, the relativistic limit of the isotropic kinetic equations (49a) and (49b) takes the simple form

ṅa k = 2π k 2 n b 1 n b 2 -2n b 1 n a k δ(k -k 1 -k 2 ) dk 1 dk 2 , ( 50a
) ṅb k = 4π k 2 n a 1 n b 2 + n a 1 n b k -n b 2 n b k δ(k 1 -k -k 2 ) dk 1 dk 2 . ( 50b
)
In Appendix C, we consider power-law spectra and show that there exist no KZ solutions with local interactions for these relativistic kinetic equations. The only power-law solutions to these equations are the relativistic limits of the thermodynamic spectra corresponding to equipartition of energy (41) or of the particle invariant (42).

We are thus left with two routes to investigate the spectral evolution of the KGMH system: In the next section we characterize the thermodynamic states reached in the long-time limit and study the wave condensation phenomenon. In Secs. VII and VIII, we come back to the nonlocal interactions that govern the out-of-equilibrium dynamics and we derive rigorously a nonlocal KZ state characterized by a nonlocal transfer of particles to a strong condensate, together with a local energy cascade in k space. This "dynamical condensation" describes the approach to the thermodynamic condensed state, and therefore makes the link between the two routes.

VI. BOSE-EINSTEIN CONDENSATION

We now study equilibration in truncated conservative systems. Nondissipative wave systems are known to exhibit a condensation phenomenon when (i) there are two conserved quantities, (ii) there is a truncation at some finite wave number K, and (iii) the spatial dimension d is large enough. In discrete systems such as the array of pendulums described in Sec. II, a truncation in k space naturally arises. In quantum systems, the Bose statistics ensures that modes with large enough energy remain unoccupied, hence an effective truncation.

Condensation of classical waves was originally reported for the Gross-Pitaevskii equation [7,12]. Further studies on condensation also considered equations of the Klein-Gordon type [13]. For all these systems, one can compute KZ solutions corresponding to an inverse cascade of the particle invariant. Condensation then naturally occurs, as the inverse cascade induces an accumulation of particles at low k. In contrast with these systems, we did not find any such KZ inverse cascade for the KGMH equation. Furthermore, the particle invariant is only an adiabatic one. In the following, we follow an approach similar to Connaughton et al. [7] to show theoretically and numerically that condensation nevertheless occurs.

We consider decaying solutions with fixed initial number of particles N and several values of the initial energy E. After some transient, we expect the system to reach the thermodynamic spectra (43), with a pair (T ,μ) corresponding to a pair (E,N ). However, we will see that there are some values of (E,N ) for which no pair (T ,μ) can be found with positive μ.

From the thermodynamic spectra, we compute the numbers N a and N b of a and b particles, as well as the energy E a and E b contained in these two fields. Assuming for simplicity that k ∈ [0,K] is a continuous variable and considering the 2D problem, we obtain

N a = 2πL 2 T K 0 k 2μ + √ 2 + k 2 dk = 2πL 2 T 2 + K 2 - √ 2 + 2μ ln √ 2 + 2μ √ 2 + K 2 + 2μ , ( 51a
)

N b = 2πL 2 T K 0 k μ + k dk = 2πL 2 T K + μ ln μ K + μ , ( 51b
)
E a = 2πL 2 T K 0 2 + k 2 k 2μ + √ 2 + k 2 dk = 2πL 2 T K 2 2 + 2μ( √ 2 -K 2 + 2) + 4μ 2 ln √ 2 + K 2 + 2μ √ 2 + 2μ , ( 51c
)
E b = 2πL 2 T K 0 k 2 μ + k dk = 2πL 2 T K 2 (K -2μ) + μ 2 ln 1 + K μ . ( 51d
)
We define the average energy per particle as

F(μ) = E N = E a + E b 2N a + N b . ( 52
)
F is independent of temperature, and can be specified at the beginning of a numerical simulation by an appropriate choice of the initial condition. However, not all values of F are compatible with the thermodynamic spectra. Indeed, we show in Fig. 7 a plot of μ as a function of F. There is a minimum value of F under which μ becomes negative: the thermodynamic spectrum is not a valid solution anymore, and wave condensation occurs. In the limit of large K, the threshold is We stress the fact that in the present 2D system, condensation occurs even in the infinite domain limit. This does not contradict the Mermin-Wagner theorem because the dispersion relation of the b waves goes to zero slow enough (ω b k ∼ k as k → 0) for the integrals (51a) and (51b) to converge at k = 0, even when μ = 0: for a given value of the energy, only a finite number of particles can be in the Rayleigh-Jeans spectrum, and the rest of the particles have to condense. By contrast, in the 2D Gross-Pitaevskii equation the dispersion relation is ω k = k 2 , and the number of particles contained in the Rayleigh-Jeans distribution involves a diverging integral for μ = 0: for a given value of the energy, an arbitrary large number of particles can be in the Rayleigh-Jeans spectrum, and condensation does not occur in 2D. For a system with dispersion relation ω k ∼ k σ in dimension d, the criterion for Bose-Einstein condensation to occur is σ < d [14].

F cr = K/3.
There is also a limiting value of F(μ) that is attained for μ → ∞. In the limit of large K, this limit is F(μ → ∞) = K/2. As one approaches this limiting value of F, the spectrum tends to the equipartition of N . Obviously, for a finite number of particles, taking the limit μ → ∞ also requires T → ∞. If the initial values of N and E are such that F > F(μ → ∞), then there is an accumulation of energy in the high k modes. This accumulation corresponds to Rayleigh-Jeans thermodynamic spectra (43) with negative values of T and μ.

To investigate numerically these different thermodynamic states, we perform numerical simulations with an initial condition φ k ∼ k α and the same initial value N 44. Large negative values of α correspond to low F and steep decreasing spectra, whereas positive α correspond to large values of F. In the following, we describe the stationary state attained after a long integration time. The size of the domain is [0,12π ] 2 and the resolution is 128 2 . After dealiasing, we keep 64 modes in each direction, i.e., we keep wave numbers k with k K = 32/6 where the step is δk = 1/6 for the two components of k. We plot the spectra of a and b modes, which we evaluate using

|a k | 2 ω (a) k | λk | 2 and |b k | 2 ω (b) k | χk | 2 .
More specifically, we show the 1D spectra S a (k) and S b (k) defined in such a way that S a;b (k)dk is the variance of the corresponding field. The four runs described in Figs. 8 to 11 illustrate, respectively, (i) a Rayleigh-Jeans state, (ii) the threshold for wave condensation, (iii) wave condensation, and (iv) accumulation at large k. In the latter simulation, the system has not settled to a steady state: the ratio F = E/N very slowly decreases in time until it reaches its critical value and a Rayleigh-Jeans state. The drift of F probably originates from the fact that N is only an adiabatic invariant: it is not conserved by the higher order nonlinearities that we discarded in the present weakly nonlinear study.

A. Nonuniform condensates

To study the condensation transition, we performed a set of numerical simulations with decreasing initial values of F: as F decreases, so does μ, until the denominator of the Rayleigh-Jeans spectrum (43) vanishes for some value of k. This happens first for a b mode since ω b k < ω a k for a given k. As a consequence, close to threshold we expect condensation of b waves, but no condensation of a waves: particles accumulate in the lowest-k b modes, which in the equivalent particle physics system corresponds to Bose-Einstein condensation of pions.

These b modes cannot have vanishing wave number because b k=0 is fixed in the leading order. Indeed, inserting the weakly nonlinear expansion ψ(x,t) = 1 + [λ(x,t) + iχ(x,t)] in the expression for the vanishing charge Q, we obtain to lowest order that d dt ( χ dx) = 0, i.e., the uniform phase mode remains constant in time, and it can be set to zero. As a consequence, in a finite periodic domain of size L in d dimensions, there are d b modes corresponding to the lowest nonvanishing frequency ω b = k = 2π/L. For such a finitesize system, condensation occurs as soon as μ = -2π/L, and it manifests as an accumulation of particles in these large-scale (but nonuniform) b modes.

As the condensate fraction increases, the weakly nonlinear assumption rapidly breaks down in the low-k modes. Indeed, the level of nonlinearity is determined by ratio of the typical intensity to the typical frequency of the waves. Hence, nonlinearity increases during condensation because b particles accumulate in the low-frequency modes. Note that conservation of the particle invariant N may even fail because of this increased nonlinearity. Nevertheless, we observe in the numerics that condensation is a robust phenomenon, that occurs even when the particle invariant is only an adiabatic one. As condensation proceeds, the enhanced nonlinearity arrests further accumulation and sharpening of the spectrum in the lowest (nonzero) frequency modes, resulting in a condensate that spreads over a few low-frequency modes. Because dynamical separation of the normal modes of the linear system is lost for strong nonlinearity, condensation of the b modes "spills" into the low-k a modes, even though the weakly nonlinear theory predicts no a-mode condensation. The resulting condensation behavior is displayed in Fig. 10.

B. Condensed fraction

To extract the condensate fraction from the numerical simulations, we first determine the temperature of the system using a fit to the high-k part of the spectra. These spectra consist of a flat isothermal level, plus a bump at low k corresponding to the condensate. By subtracting the isothermal part of the spectra before integrating over k we access the number of condensed particles. The condensed fraction is the ratio of the number of condensed particles to the total number of particles. We plot it in Fig. 12. From the numerics, the threshold for condensation is estimated to be F 2.1. We can predict this value accurately if we take into account the discreteness of the lattice of wave numbers. Replacing the integrals in Eqs. (51a) to (51d) by a sum over the discrete lattice of wave numbers, we obtain a critical value for condensation F cr = 2.06, hence, a very good agreement between the theory and the numerics.

The condensed fraction can be predicted theoretically as follows: Under the threshold for condensation, the chemical potential vanishes and the thermodynamic spectrum of b particles diverges for small k; slightly below the onset F cr , a macroscopic number N (cond) b of b particles condenses, while the a particles and the remaining b particles are (approximately) in energy equipartition, i.e., they have μ 0 thermodynamic spectra. The noncondensed particles carry (almost) all the energy E of the system, and their number is

2N a + N b -N (cond) b
. The ratio of these two quantities is F cr , from which we deduce the condensed fraction of b particles

N (cond) b N = F cr -F F cr . ( 53
)
This prediction is in good agreement with the numerical simulations (see Fig. 12).

VII. STATISTICAL DESCRIPTION OF THE DECAY INSTABILITY

Kinetic equations allow to compute the out-of-equilibrium evolution of some initial distribution of waves. In particle physics, such an initial distribution can be given in some cosmological model, or it can correspond to particles appearing at t = 0 in an experiment. As we mentioned earlier, for the KGMH system the transfers governed by the kinetic equations can be very nonlocal in k space. As an example of such nonlocal transfers, we come back to the decay instability described in Sec. IV, which we now study for isotropic ensembles of waves.

A. Decay instability for an isotropic ensemble of a waves

Consider a 3D isotropic ensemble of a waves with wave number κ. These waves spontaneously induce some b waves through the decay instability. A key question is to know the distribution of b waves that is produced by this instability: What is the wave number (or energy) of the b waves? In the particle physics analogous system, we use the kinetic equation to predict the distribution of pions that arises from the decay of σ mesons with wave number κ.

The initial distribution is

n a k = N a 4πκ 2 L 3 δ(k -κ), ( 54a
)
n b k = 0. ( 54b
)
As already mentioned in Sec. V C, this distribution is a steady solution to the kinetic equations because n b k = 0. However, it is an unstable solution, as we can see by considering an infinitesimal perturbation n b k (t) 1. Inserting (61a) into kinetic equation (49b) and retaining only the linear terms in the infinitesimal perturbation n b k , we obtain

ṅb k = N a k 2 L 3 κ √ 2 + κ 2 n b k + n b √ 2+κ 2 -k if k ∈ I ≡ √ 2 + κ 2 -κ 2 ; √ 2 + κ 2 + κ 2 , (55a)
ṅb k = 0 otherwise, (55b)
where the condition in (55a) arises from the integration domain D b (k). We write the same equation at wave number √ 2 + κ 2k to obtain the following 2 × 2 system of linear ODEs:

ṅb k = N a k 2 L 3 κ √ 2 + κ 2 n b k + n b √ 2+κ 2 -k , ( 56a
) ṅb √ 2+κ 2 -k = N a L 3 ( √ 2 + κ 2 -k) 2 κ √ 2 + κ 2 n b k + n b √ 2+κ 2 -k , (56b) if k ∈ I, and ṅb k = ṅb √ 2+κ 2 -k = 0 otherwise.
Consider an eigenmode (n b k ,n b √ 2+κ 2 -k ) = ( nb k , nb √ 2+κ 2 -k )e σ t
, where σ is the growth rate. Substitution into the system of equations yields a 2 × 2 system in nb k and nb √ 2+κ 2 -k , that has nontrivial solutions when the determinant vanishes. This gives the following expression for the reduced growth rate σ = σ κ

√ 2 + κ 2 L 3 /N a : σ (k) = 1 ( √ 2 + κ 2 -k) 2 + 1 k 2 if k ∈ I , (57a) σ (k) = 0 otherwise. ( 57b
)
Plots of σ (k) are provided in Fig. 13, in the limiting cases κ 1 and κ 1. For κ 1, I is a narrow interval centered around

√ 2+κ 2 2 1 √
2 , and b waves appear at this wave number: there is a strong nonlocal transfer of massive particles with κ 1 into massless particles with k

1 √ 2 .
For κ 1, the growth rate has a maximum at k κ and a maximum at

k = √ 2 + κ 2 -κ 1 κ
1: massless particles will appear preferentially at these two wave numbers, hence, a strong nonlocal transfer from a particles with k = κ into b particles with k = 1 κ 1. These predictions for isotropic initial distributions of a waves are qualitatively compatible with the numerical results of Sec. IV on the decay of a single a mode. If the initial a mode has wave number κ 1, b waves appear at both wave numbers k 1 and k κ (see Fig. 3), whereas if the initial a mode has wave number κ 1, then b waves appear at wave number k 1 √ 2 (see Fig. 4).

B. Reaction yield of the decay instability

Starting from an isotropic distribution of a particles only with wave number κ, together with weak background noise on the a and b fields, the decay instability sets in and b waves appear through nonlocal transfers. The intensity of the a waves decreases accordingly, to ensure conservation of the energy and particle invariants. When b waves become strong enough, the quadratic term in ñb k of kinetic equation (49a) cannot be neglected anymore. This term modifies the spectral distribution of a waves, the system enters the nonlinear regime, and eventually it reaches thermodynamic equilibrium, with both spectra of a and b waves spreading over the entire range of wave numbers.

We can thus determine the saturated state of the decay instability using the thermodynamics described in Sec. VI. Indeed, the initial distribution (61a) corresponds to an energy per particle

F = √ 2 + κ 2 2 . ( 58
)
An a wave can be seen as a molecule made of two b atoms, and the decay instability corresponds to the reaction a → b + b.

In the saturated state, one can compute the yield η of this reaction using thermodynamics: Out of the initial N number of b atoms (two per a wave, and one per b wave), what fraction η decays into free b particles? The reaction yield is given by the equilibrium value of N b /N , and it depends on the initial wave number κ of the a waves. We now compute it in 2D, as a function of κ.

For F > F cr , no condensation occurs, and the yield is simply given by

η = N b (μ) 2N a (μ) + N b (μ) = K + μ ln μ K+μ K + μ ln μ K+μ + 2 √ 2 + K 2 - √ 2 + 2μ ln √ 2+2μ √ 2+K 2 +2μ , ( 59
)
where the chemical potential μ is related to κ through the relation (58), i.e., κ = 4F(μ) 2 -2 where F(μ) is given by (52).

For F < F cr , condensation of b waves occurs and the spectra consist of thermodynamic spectra with μ = 0, together with b particles condensed at low k. The reaction yield becomes

η = N b (μ = 0) N + N (cond) b N = F N b (μ = 0) E a (μ = 0) + E b (μ = 0) + F cr -F F cr = F K + F cr -F F cr = 1 + F 1 K - 1 F cr = 1 + 2 + κ 2 1 K - 1 F cr . ( 60
)
In Fig. 14, we plot the yield given by expressions (59) and (60) as a function of the initial wave number κ of the a particles, for a maximum wave number K 1. In this limit, F cr K/3, and we observe the following behavior: the reaction yield is a 1 (here K = 100). In the decreasing part, some fraction of the produced b particles condenses. In the increasing part, such condensation does not occur. nonmonotonous function of κ. It is close to unity for κ 1, i.e., the reaction is almost complete, with most of the a particles eventually decaying into b particles. In this regime, a large fraction of the produced b particles form a condensate. As κ increases the yield decreases, with a minimum attained at the threshold of wave condensation κ cr = √ 4F 2 cr -2 2F cr 2K/3. For κ > κ cr , the yield increases with κ. For κ = K, the yield reaches 0.5, i.e., half of the b "atoms" only are free b waves, the other half remaining in the form of a waves.

VIII. DYNAMICAL CONDENSATION AS A NONLOCAL KOLMOGOROV-ZAKHAROV CASCADE

In this section, we describe the dynamical process by which condensation is achieved. We first show that ensembles of b waves at low-enough wave numbers are steady solutions to the kinetic equations. We refer to such solutions as "condensates." We consider an infinitesimal perturbation of the aand b-wave spectra and study the linear stability of a condensate at k = 0. This analysis indicates that there is a rapid transfer of particles from the perturbation to the condensate. However, this transfer is only partial, and some marginally stable perturbations remain. We study the evolution of these remaining perturbations through a weakly nonlinear expansion and derive a kinetic equation that governs the nonlinear evolution of the perturbation to the condensate. This kinetic equation admits a KZ solution, with a net energy flux towards high wave numbers. In an infinite system, this solution corresponds to a condensate that increases in amplitude, while the energy of the remaining waves cascades to larger and larger wave numbers. In a finite system, the cutoff in wave number arrests this cascade and the system eventually reaches thermodynamic equilibrium.

A. Steady solutions with b waves only

Consider the kinetic equations for n a k = 0, and n b k = 0, and further restrict attention to distributions of b waves at low wave numbers: we assume that there is a maximum wave number K that supports the b waves, i.e., n b K = 0, and n b k = 0 for k > K. The evolution of this initial distribution of waves is governed by the kinetic equations (49a) and (49b). At t = 0 + , only the quadratic terms in n b k may be nonzero. In the righthand side of Eq. (49a), the quadratic term involves triads of wave numbers

[k 1 ,k 2 ,k = (k 1 + k 2 ) 2 -2 ] such that n b k 1 = 0, n b k 2 = 0, k 1 k 2 > 1/2, k 1 + k 2 √ 2.
Such k 1 and k 2 cannot be found if K 1 √ 2 ; provided the distribution of b waves is restricted to wave numbers smaller than 1 √ 2 , the right-hand side of (49a) vanishes. We assume that this condition is satisfied by the initial distribution, and we consider the right-hand side of kinetic equation (49b). Because n a k = 0, only the quadratic term in n b k may be nonzero. However, it is nonzero only if we can find two wave numbers k and k 2 such that n b k = 0, n b k 2 = 0, and kk 2 1/2. But, the latter condition cannot be satisfied if K < 1 √ 2 , hence, the right-hand side of (49b) vanishes.

To summarize, we have shown that an initial distribution of b waves only (n a k = 0) is a steady solution to the kinetic equations (49a) and (49b) provided these waves have only wave numbers smaller than 1 √ 2 .

B. Dynamics of waves on a condensate

For simplicity, we consider a b-wave condensate at very low wave number, which we simply write as A δ(k), and we denote as ña (k) and ñb (k) the remainders of the wave spectra. Hence,

n a k = ña k , ( 61a
)
n b k = Aδ(k) + ñb k . ( 61b
)
Inserting this decomposition into the 3D kinetic equations for k = 0 leads to

ṅa k = 4πA k √ 2 + k 2 ñb √ 2+k 2 -ña k + 2π kω a k D a (k) δ ω abb k12 ñb 1 ñb 2 -2 ñb 1 ña k dk 1 dk 2 , (62a) ṅb k = 4πA k 2 ña √ k 2 -2 -ñb k + 4π k 2 D b (k) k 1 ω a 1 δ ω abb 12k ña 1 ñb 2 + ña 1 ñb k -ñb 2 ñb k dk 1 dk 2 . ( 62b
)
The terms proportional to A arise from the interaction between the condensate and the remainders, while the integral terms correspond to interactions of the remainders. The latter terms are exactly the collision integrals of the original kinetic equations, with n replaced by ñ. We stress the fact that if ña k and ñb k are the isothermal distributions (41), we obtain a steady solution to (62a) and (62b), i.e., both the terms proportional to A and the integral terms vanish: a condensate at k = 0 together with remainders in isothermal equilibrium constitute a steady solution to the kinetic equations, in perfect agreement with the thermodynamic treatment of Bose-Einstein condensation.

To study the dynamical process of condensation, we perform a multiple-scale expansion of Eqs. (62a) and (62b) to describe the evolution of remainders much weaker than the condensate ñ A. Consider a small parameter 1 and A = O(1), and expand the remainders as ña k = ña;0 (k,t,T ) + 2 ña;1 (k,t,T ) + . . . , (63a) ñb k = ñb;0 (k,t,T ) + 2 ñb;1 (k,t,T ) + . . . ,

where the slow time is T = t.

Inserting this decomposition into (62a) and (62b) gives at order O( ) the linear stability analysis for the condensate

∂ t ña;0 k = 4πA k √ 2 + k 2 ñb;0 √ 2+k 2 -ña;0 k , ( 64a
)
∂ t ñb;0 k = 4πA k 2 ña;0 √ k 2 -2 -ñb;0 k . ( 64b
)
Considering the first of these two equations at wave number √ k 2 -2 instead of k leads to the following 2 × 2 system of first-order linear ODEs:

∂ t ña;0 √ k 2 -2 = 4πA √ k 2 -2 k ñb;0 k -ña;0 √ k 2 -2 , ( 65a
)
∂ t ñb;0 k = 4πA k 2 ña;0 √ k 2 -2 -ñb;0 k . ( 65b
)
Solutions to this system of ODEs can be sought in the form ( ña;0 √ k 2 -2 , ñb;0 k ) ∼ e σ t , which leads to the following growth rates:

σ 1 = 0 and σ 2 (k) = -4πA 1 k 2 + 1 k √ k 2 -2 < 0. ( 66
)
Decomposing on the corresponding eigenmodes, we obtain

ña;0 √ k 2 -2 = C 1 (k,T ) + C 2 (k,T )e σ 2 (k)t , ( 67a
) ñb;0 k = C 1 (k,T ) - √ k 2 -2 k C 2 (k,T )e σ 2 (k)t . (67b)
Because σ 2 < 0, the contribution in C 2 rapidly decays on the fast time t, and after a short transient,

ña;0 √ k 2 -2 = ñb;0 k ≡ n k (T ). ( 68
)
One can check that the C 2 part of the spectra has vanishing total energy at any time. Its rapid decay occurs at constant energy, as it should. By contrast, this C 2 contribution carries a nonzero number of particles

N 2 = C 2 (k) √ k 2 -2ke σ 2 (k)t dk.
The decay of the C 2 contribution on the fast time t therefore corresponds to a rapid transfer of particles between the remainders and the condensate. These rapid transfers of particles ensure that (68) is satisfied adiabatically. The dynamics of the system is then characterized by a single slowly evolving spectrum n k (T ) from which both spectra of a and b particles can be deduced. We pursue the expansion to determine the kinetic equation that governs the evolution of n k (T ).

To order O( 2 ), we consider a linear combination of Eq. (62a), considered at wave number √ k 2 -2, and of Eq. (62b) to obtain

∂ t k 2 -2 ña;1 √ k 2 -2 + k ñb;1 k = -(k + k 2 -2) ∂ T n k + 2π k k- √ k 2 -2 2 k 1 k+ √ k 2 -2 2 n k 1 n k-k 1 -2n k 1 n k dk 1 + 2 k 1 k+ 1 2k n k 1 n k 1 -k + n k 1 n k -n k 1 -k n k dk 1 . ( 69
)
To obtain this form for the second integral, we have made the change of variables

√ k 2 1 + 2 → k 1 in the integrand.
The solvability condition demands that the right-hand side vanish, which leads to the following kinetic equation:

∂ T n k = 2π k(k + √ k 2 -2) × k- √ k 2 -2 2 k 1 k+ √ k 2 -2 2 n k 1 n k-k 1 -2n k 1 n k dk 1 + 2 k 1 k+ 1 2k n k 1 n k 1 -k + n k 1 n k -n k 1 -k n k dk 1 . ( 70
)
In the relativistic limit k 1, the equation becomes scale invariant:

∂ T n k = π k 2 k 1 k n k 1 n k-k 1 -2n k 1 n k dk 1 + 2 k 1 k n k 1 n k 1 -k + n k 1 n k -n k 1 -k n k dk 1 . ( 71
)
This equation conserves the total energy, but it does not conserve the number of particles in the remainders. Indeed, these remainders rapidly exchange particles with the condensate, and only the total number of particles contained in both the remainders and the condensate is conserved. In Appendix C 2, we determine the steady power-law solutions to this equation n k ∼ k x . As already mentioned, one such solution to the kinetic equation is the isothermal spectrum n k ∼ k -1 that coexists with the condensate when wave condensation occurs. But, this kinetic equation also admits a solution of the KZ type, with exponent x = - 3 2 . This turbulent spectrum carries a flux P of energy through the scales of the system. It reads as (see Appendix C 2 for details)

n k = P 1 2 k -3 2 4π √ 2π -8 ln 2 . ( 72
)
In this cascading state, the flux of energy is due to local interactions in k space, while nonlocal interactions continuously transfer particles from the remainders to the condensate to satisfy the adiabatic constraint n k = n a k = n b k [one could highlight the corresponding slow growth of the condensate by considering the O( 2 ) equations at k = 0 and introducing an evolution of the amplitude A of the condensate on an even slower time scale 2 t].

IX. DISCUSSION

We have developed a weak-turbulence theory for the Klein-Gordon equation with a Mexican-hat potential (1), a universal nonlinear equation which displays spontaneous symmetry breaking. The fundamental state of this equation is a uniform field at the minimum of the Mexican-hat potential. The weakturbulence theory describes the interaction of random weakly nonlinear waves propagating over this uniform background state. Such wave disturbances are of two kinds, with two branches in the dispersion relation: in the particle physics context, these branches correspond to a massive and a massless particle field, and the KGMH equation provides the simplest particle physics model for a system of interacting σ mesons and pions. In the weak-nonlinearity limit, the dominant interaction between these fields is a three-wave process by which two massless waves can merge into a massive one or, inversely, one massive wave can decay into two massless ones. This process conserves energy and an equivalent of the number of particles.

We have derived the wave-turbulence kinetic equations describing the coupled evolution of the spectra of massive and massless waves under this three-wave interaction. They admit a two-parameter family of thermodynamic equilibrium spectra of the Rayleigh-Jeans type. When the initial condition corresponds to little energy per particle, the thermodynamic state reached in the long-time limit displays massless-wave condensation, the classical analog of Bose-Einstein condensation.

The kinetic equations admit no Kolmogorov-Zakharov (KZ) cascading states in the usual sense. Instead, the out-ofequilibrium dynamics is dominated by nonlocal interactions in scale space. We have highlighted two cases of such nonlocal evolution. The first one is when a massive mode spontaneously decays into massless waves, and the second one is when small-scale a and b fields evolve in presence of an intense large-scale b condensate. In the latter regime, some particles are rapidly transferred to the condensate, and the small-scale remaining waves follow a single reduced kinetic equation which admits a KZ solution corresponding to an energy cascade.

These different regimes are illustrated in the 3D numerical simulation reported in Fig. 15

S λ (k) = ( L 2π ) 3 4πk 2 |λ k | 2 and S χ (k) = ( L 2π ) 3 4πk 2 |χ k | 2 .
The simulation starts off with a distribution of massive waves only, together with a weak level of massless waves. The initial spectrum S λ is flat up to wave number k = 6, i.e., the massive waves are mostly relativistic [see Fig. 15(a)]. The initial evolution of this spectrum corresponds to a fast decay instability. As discussed in Sec. VII, the decay of relativistic massive particles with wave number k a produces massless waves with k b k a and with k b 1. The latter waves accumulate in a strong condensate of massless particles. The subsequent evolution corresponds to waves superposed to a strong condensate, the dynamics of which we presented in Sec. VIII: the spectra of the small-scale remaining a and b waves rapidly become identical [see Fig. 15(b)], and follow the reduced kinetic equation (71). After some transient, we obtain the quasistatic spectral distribution plotted in Fig. 15(c), which corresponds to the KZ energy cascade (72). Indeed, we computed the energy leaving the domain k ∈ [1; 10] per unit time to determine the value of the energy flux. We obtained P 7 × 10 -4 . In the relativistic regime, the KZ spectra corresponding to (72) are given by kS λ (k) = kS χ (k) = 4πk 2 n k , i.e., S λ (k) = S χ (k) = P

1 2 k -1 2 √ 2π -8 ln 2 . ( 73
)
This theoretical prediction is plotted in Fig. 15(c). It has no fitting parameters and is in good quantitative agreement with the spectra obtained numerically in the inertial range. This confirms that the system displays a nonlocal Kolmogorov-Zakharov cascading state, where nonlocal transfers of particles with the condensate coexist with a local energy cascade. This state is therefore the intermediate asymptotics that describes the evolution between the initial distribution of waves and the thermodynamic equilibrium attained in the long-time limit. The convergence of these integrals depends on the value of x, C a , and C b . For x > 1 2 , the right-hand side of (C4) diverges unless C a = C b . For x ∈]0; 1/2] and x ∈ [-1/2; 0[, the righthand side of (C4) diverges. For x = 0 it converges only if C b = 2C a , for which it vanishes. The right-hand side of (C3) also vanishes for (x = 0, C b = 2C a ). The corresponding solution is the equipartition of particle number N . For x ∈] -2; -1/2[, the right-hand side of (C4) converges provided C a = C b . When this condition is met, the right-hand side of (C3) also converges, and the equations become

ṅa k = 2πC 2 k 2x-1 - 2 1 + x + 2 (1 + x) (2 + 2x) , ( C5a
)
ṅb k = 2πC 2 k 2x-1 - 2 1 + x + 2 (-1 -2x) (1 + x) (-x) , ( C5b
)
where C = C a = C b . We plot in Fig. 16 the arguments of the square brackets in the right-hand side of these two equations, i.e., we plot, respectively, The first one indicates if the spectrum of energy is stationary, whereas the second one indicates whether the spectrum of particle number N is stationary. We see that ṅa k + ṅb k vanishes at x = -1 but also at x = -3 2 . On dimensional grounds (see e.g. [START_REF] Nazarenko | Wave Turbulence[END_REF]), one can see that the latter value corresponds to the KZ spectra for an energy cascade. However, this spectrum is not a stationary solution to the full system of equations because ṅa k and ṅb k do not vanish independently at x = -3 2 .

The combination 2 ṅa k + ṅb k has a double zero at x = -1. Indeed, this index corresponds both to energy equipartition and, on dimensional grounds, to a constant flux of N . Still, it remains unclear whether it is possible to observe a stationary state with such a nonzero flux of N : the value x = -1 is a single zero of ṅa k + ṅb k , which corresponds to energy equipartition. It is thus possible that ṅa k + ṅb k is not zero when there is a flux of N .

We want to solve this equation in a periodic domain [0,2π ] d . The parameter g > 0 allows to put the crossover wave number anywhere in k space, i.e., the dispersion relation for a waves is now ω a k = 2g + k 2 . For better accuracy, one wants to integrate exactly the rapidly oscillating part of the solution. Here, we write the equation as ψ tt = ψ + NL with NL = (|ψ| 2g)ψ .

(D2)

If we drop the term NL, a solution in Fourier space is ψ k = A + e ikt + A -e -ikt . Let us use the method of variation of constants: we write ψ k = A + (t)e ikt + A -(t)e -ikt and impose the constraint A + (t)e ikt + A -(t)e -ikt = 0. Inserting this decomposition into Eq. (D2) we obtain

A + = i 2k e -ikt NL k , ( D3a
)
A -= - i 2k e ikt NL k , ( D3b
)
where NL k is the Fourier amplitude of (|ψ| 2g)ψ for wave vector k. We integrate Eqs. (D3a) and (D3b) using an Adams-Bashforth scheme with time step dt:

A + (
For nonzero k, the Fourier amplitude of the field is then ψ k (t) = ψ + k (t) + ψ - k (t). The mode k = 0 must be taken care of separately. Using a standard finite-difference approximation to the second order time derivative, we get ψ 0 (t + dt) = 2ψ 0 (t)ψ 0 (tdt)dt 2 NL 0 .

(D6)

We compute the nonlinear term using a standard pseudospectral method with dealiasing following the 1/2 rule. A time step consists of the following: the nonlinear term is computed in real space. The field ψ is computed as the inverse Fourier transform of ψ k and the nonlinear term (|ψ| 2g)ψ is transformed back to Fourier space to get NL k (t). Then, Eqs. (D5a), (D5b), and (D6) are iterated to obtain ψ k at time t + dt.

FIG. 2 .

 2 FIG. 2. (Color online) Two types of waves propagate in this system: branch a corresponds to radial oscillations of the pendulums, whereas branch b corresponds to azimuthal oscillations of the pendulums, along the minimum of the Mexican-hat potential.

FIG. 3 .

 3 FIG. 3. (Color online) Destabilization of a small-scale massive wave into massless ones. The left-hand panels are the λ fields, the middle panels are the χ fields, and the right-hand panels are the 1D spectra of λ (blue line) and χ (red dashed line) as a function of k.Here, L = 20π , and from top to bottom: t = 0, t = 30, t = 60, and t = 120. The unstable triads involve a small-scale massless mode and a large-scale one.

FIG. 4 .

 4 FIG. 4. (Color online) Destabilization of a large-scale massive wave into massless ones. The left-hand panels are the λ fields, the middle panels are the χ fields, and the right-hand panels are the 1D spectra of λ (blue line) and χ (red dashed line) as a function of k. Massless waves appear at the crossover scale k 1. Here, L = 20π , and from top to bottom: t = 0, t = 120, t = 270, and t = 570.

D 1 √ ω a k ω b 2 ω b 3

 13 . Isotropic systems Consider the class of isotropic spectra n a k ≡ n a k and n b k ≡ n b k , where k = |k|, and integrate out the angular variables describing the directions of k 1 and k 2 in the kinetic equations (38a) and (38b). Because V k 12 ≡

FIG. 5 .

 5 FIG. 5. Integration domain D a (k) for k < 1/ √ 2 (left-hand panel) and k > 1/ √ 2 (right-hand panel).

)FIG. 6 .

 6 FIG. 6. Integration domain D b (k) for k < 1/ √ 2 (left-hand panel) and k > 1/ √ 2 (right-hand panel).

FIG. 7 .

 7 FIG. 7. (Color online) The chemical potential μ as a function of the "energy per particle" F = E/N . The usual Rayleigh-Jeans spectrum is a solution for intermediate values of F only. For low F condensation occurs at low k, whereas for large F energy and particles accumulate at the cutoff wave number K.

  FIG. 8. (Color online) α = -0.5: in the stationary state, F = 2.42 ∈ (K/3,K/2), so we expect a Rayleigh-Jeans distribution. Top: spectra of a and b modes (left), and spectra normalized by a Rayleigh-Jeans spectra with μ = 1.67 and T = 7 × 10 -4 (right). Bottom: snapshots at final time of the λ and χ fields.

  FIG. 9. (Color online) α = -0.8: in the stationary state, F = 2.12, which is close to the critical value (approximately K/3) for wave condensation on this K = 32/6 grid. Top: spectra of a and b waves (left), and spectra normalized by pure energy equipartition spectra with T = 3.5 × 10 -4 (right). Bottom: snapshots at final time of the λ and χ fields.

  FIG. 10. (Color online) α = -1.3: in the stationary state, F = 1.12, which is below the value for wave condensation. Top: spectra of a and b waves (left), and spectra normalized by energy equipartition spectra with T = 3.3 × 10 -5 . The flat part at high k is in thermal equilibrium, and the bump at low k indicates condensation. Bottom: snapshots at final time of the λ and χ fields.

FIG. 11

 11 FIG. 11. (Color online) α = -2: after 6000 time units, F = 3.08, which is above the critical value for high-k accumulation. Top left: spectra of a and b waves. Top center: spectra normalized by equipartition of N , i.e., Rayleigh-Jeans spectra in the limit (T ,μ) → (∞,∞), with T /μ = 1.5 × 10 -4 . The flat part at low k is in equipartition of N , and the bump at high k indicates accumulation. Top right: spectra normalized by thermodynamic spectra with negative temperature T = -1.3 × 10 -3 and chemical potential μ = -10. Bottom: snapshots after integration over 7200 time units of the λ and χ fields.

  FIG. 12. (Color online) Condensate fraction as a function of the energy per particle E/N . We show the condensate fraction of a particles (•), b particles ( ), and of 2N a + N b ( ). The dashed line is the theoretical prediction (53).

FIG. 13 .

 13 FIG. 13. (Color online) Reduced growth rate of b waves as a function of wave number k, for a waves at wave number κ. The dashed curve is ( √ 2 + κ 2k) -2 + k -2 . Left-hand panel: for κ 1, b waves grow at wave number k 1/ √ 2. Right-hand panel: for κ 1, the growth rate has two strong maxima, one at k 1/κ and one at k κ.

FIG. 14

 14 FIG. 14. (Color online)Reaction yield η for d = 2 as a function of the initial wave number κ of the a waves, for a maximum allowed wave number K 1 (here K = 100). In the decreasing part, some fraction of the produced b particles condenses. In the increasing part, such condensation does not occur.

√ 6 × 6 20.

 66 : the domain is [0,L] 3 with L = 2π √ 6 15.4, with a resolution of 192 3 . The wave vector k is on the 3D grid 1 [-96,96] 3 , but after dealiasing we keep only the modes with k 48 √ We plot the 1D wave spectra

FIG. 15

 15 FIG. 15. (Color online) Out-of-equilibrium spectra S λ (blue solid line) and S χ (red dashed line) in a 3D numerical simulation of the KGMH equation. Panel (a) shows the initial condition [S χ (k) 2 10 -3 is out of the figure]. Panel (b) corresponds to t = 40: after a fast transient, part of the massive waves decayed into massless waves. A strong condensate of b waves appears, and the small-scale remaining spectra follow the adiabatic constraint S λ S χ . Panel (c) corresponds to t = 475: energy cascades from low to high wave numbers, in quantitative agreement with the KZ prediction (73), shown as a thick solid line.

FIG. 16

 16 FIG. 16. (Color online) Collision integrals as a function of the index x. (a)

ṅa k 2πC 2 k 2x-1 and ṅb k 2πC 2 k 2x- 1 .

 1 Both have a single zero at x = -1 for the range x ∈] -2; -1/2[: this can be either the equipartition of wave energy, or the KZ spectra for a cascade of N .It is insightful to also plot the linear combinations ṅa

  t + dt) = A + (t) + i dt 4k e -ikt [3NL k (t)e ikdt NL

	i dt 4k [3NL i dt 4k [-3NL

k (tdt)], (D4a)

A -(t + dt) = A -(t) + i dt 4k e ikt [-3NL k (t) + e -ikdt NL k (tdt)]. (D4b) Let us define ψ + k (t) = A + (t)e ikt and ψ - k (t) = A -(t)e -ikt : ψ + k (t + dt) = e ikdt ψ + k (t) + k (t)e ikdt NL k (tdt)] ,

(D5a)

ψ - k (t + dt) = e -ikdt ψ - k (t) + k (t) + e -ikdt NL k (tdt)] .

This definition of T is common in wave turbulence. In the energy equipartition state, it corresponds to an energy k B T per oscillatory degree of freedom, where k B is a dimensionless Boltzmann constant k B = (2π ) d in d dimensions.

ACKNOWLEDGMENTS

We thank J.-P. Blaizot for insightful discussions. During the time this work was conducted, S. Nazarenko was supported by Labex PALM ANR-10-LABX-0039 at SPEC laboratory, CEA, France.

APPENDIX A: THREE-WAVE RESONANCES

The interaction Hamiltonian (20) corresponds to threewave processes of type a + a + a → 0, a + a → a, a + b + b → 0, a + b → b, and b + b → a. Out of these processes, the only ones allowed are those that can satisfy both the wave number and the frequency resonance conditions. Clearly, these processes cannot be resonant because the sum of three positive frequencies cannot be zero. Thus, we neglect the corresponding terms in the Hamiltonian.

a + a → a process

The frequency resonance condition for this process is ω a k 1 + ω a k 2 = ω a k 3 . After the substitution 2 . Taking the square of both sides, we obtain

. This equation has no solutions because the left-hand side is strictly greater than k 1 k 2 whereas the right-hand side is less than k 1 k 2 . Thus, we neglect the corresponding terms in the Hamiltonian.

a + b → b process

For this process, the resonance conditions give

Squaring both sides yields

Again, this equation has no solution because the left-hand side is strictly greater than k 1 k 2 whereas the right-hand side is less than k 1 k 2 , so we neglect the corresponding terms in the Hamiltonian.

b + b → a process

For this process, the resonance conditions give

Squaring both sides, we obtain

where θ is the angle between k 1 and k 2 . We have

i.e., such if

This condition predicts nonlocal interactions if the initial spectrum of a waves is at k 1; then, the initial evolution generates b waves at k ≈ 1/ √ 2. Indeed, k 1 implies k 1 ≈ k 2 and cos θ ≈ -1, which according to (A1) gives

Conversely, two short b waves can produce a long a wave through a wave-beating process. On the other hand, if the a particles are relativistic then the condition (A2) is not restrictive: for k,k 1 ,k 2 1 we have cos θ ≈ 1, which is typical of acousticlike systems; the wave vectors are almost collinear in all resonant triads.

APPENDIX B: WAVE-TURBULENCE DERIVATION

Wave-turbulence deals with statistical ensembles of waves. For such wave fields, a small amount of dispersion rapidly leads to the complex amplitudes of the different waves being uncorrelated. We can therefore consider the amplitudes and the phases of the waves to be random independent variables to start with. Such a random phase and amplitude approximation is a useful shortcut to obtain the wave kinetic equations. The brackets . . . denote an ensemble average over these random phases and amplitudes.

Separating the time scales

Consider the interaction representation variables

where the small parameter 1 has been introduced to study the weakly nonlinear behavior. These variables are steady in the absence of nonlinearity, and they evolve slowly for weak nonlinearity, with ċk = -i

where where the inequality must be satisfied for both ω = ω a and ω b .

Seek the solutions at time T via a regular expansion in :

Following, we neglect the orders of higher than two. Substitution into the evolution equations and (B1b) gives at leading k = ḋ(0) k = 0, i.e.,

At the next order we have

where

Doing one more step of the recursion, we get

Substituting here c (1) and d (1) from (B4a) and (B4b) and integrating over time yields

where

Introduce the wave intensities

and expand them to second order in ,

Phase averaging

We consider the phases of the waves as random variables that are uncorrelated for different wave vectors. When phase averaging the expressions (B8a) and (B8b), the linear order in drops out, as usual in wave turbulence. Denoting the average over those phases as . . . φ , we obtain c (1)

where

We can further average over random initial amplitudes: the intensities at different wave numbers are independent random variables, for instance,

L → ∞ and T → ∞ limits: Kinetic equations

In wave turbulence, the limit L → ∞ is taken before the limit T → ∞, the latter being equivalent to the → 0 limit. For the limit L → ∞, we follow the standard rules to switch to a continuous description

and replace the Kroenecker deltas with Dirac deltas:

For the T → ∞ limit, we take into account the following limits:

and

Introduce the wave spectra

and performing the infinite-domain and long-time limits, we arrive at the kinetic equations

Recall that here δ k 12 , etc., denote Dirac deltas rather than Kroenecker deltas.

APPENDIX C: ISOTROPIC POWER-LAW SPECTRA IN THE RELATIVISTIC LIMIT 1. Calculation of the collision integral for power-law spectra

We consider the case d = 3, and we assume power-law spectra of the form

At first sight, we may be tempted to consider a more general form in which the index is different for n a k and n b k , but that proves unfruitful: different indices would mean different degrees of homogeneity of the different contributions to the right-hand side of the kinetic equations and they could not sum up to zero. Second, to find the steady-state power-law solutions, one might like to use the so-called Zakharov transformation, as is customary in WT theory [START_REF] Falkovich | Kolmogorov Spectra of Turbulence[END_REF][START_REF] Nazarenko | Wave Turbulence[END_REF]. However, we will see in the following that in our system all integrals can be evaluated analytically for any exponent x, and, therefore, the Zakharov transformation is unnecessary.

Inserting (C1) in Eq. (50a) and integrating over k 2 , one gets

Changing variable to u = k 1 /k, we obtain

Performing the same steps to Eq. (50b) yields

KZ cascade on a condensate

The remainders that evolve together with a strong condensate follow the kinetic equation (71) in the relativistic limit, where the a and b waves have the same spectrum n k (T ). Let us look for power-law solutions n k = Ck x to this equation. Substituting into (71) and changing the integration variable to u = k 1 /k leads to

where

The function J (x) vanishes for x = -1 and -3 2 , the former corresponding to energy equipartition and the latter being a KZ-type solution. Let us compute the energy flux P in this KZ state: the kinetic equation is a conservation equation for the total energy, that we can write in terms of the energy spectrum (ω a + ω b )4πk 2 n k = 8πk 3 n k . Considering powerlaw solutions n k = Ck x , ∂ T (8πk 3 n k ) = -

The flux P of energy going from k < k 0 to k > k 0 per unit time is obtained by integrating the right-hand side expression from k = 0 to k 0 :

For the KZ exponent x = -3 2 , this flux is independent of k 0 , indicating a constant flux of energy cascading through the scales of the system. Indeed, both the numerator and the denominator of (C9) vanish at x = - 3 2 , and Taylor expansion of J (x) to first order around this value leads to P = 32π 2 C 2 (π -4 ln 2). We can therefore express C in terms of P, which leads to the KZ solution