Verification of Wave Turbulence Theory in the Kinetic Limit

Using the 1D Majda-McLaughlin-Tabak model as an example, we develop numerical experiments to study the validity of the Wave Kinetic Equation (WKE) at the kinetic limit (i.e., small nonlinearity and large domain). We show that the dynamics converges to the WKE prediction, in terms of the closure model and energy flux, when the kinetic limit is approached. When the kinetic limit is combined with a process of widening the inertial range, the theoretical Kolmogorov constant can be recovered numerically to a very high precision.


I. INTRODUCTION
Wave Turbulence (WT) describes the out-ofequilibrium statistical dynamics of multi-scale wave systems.The theory of WT has been successfully applied to various physical contexts, including ocean surface waves [1,2], internal gravity waves [3,4], quantum turbulence [5], and gravitational waves in the early universe [6].For a given system, a statistical closure model can be developed that connects the high-order correlators of the wave field to pair correlators.When the closure model is taken in the kinetic limit, i.e., infinitesimal wave amplitude in an infinite domain, a wave kinetic equation (WKE) can be derived which describes the spectral evolution via a Boltzmann-like collision integral over wave-wave interactions.
One of the most important features of the WKE is that it yields stationary power-law solutions with constant flux, known as Kolmogorov-Zakharov (KZ) spectra.For direct cascades, the general form of the KZ solution of wave action spectrum can be written as n k = CP α k γ , where k is the wavenumber, P is the energy flux, and α = 1/2 and α = 1/3 for systems with three and fourwave resonances, respectively.C and γ are constants that can be calculated as a part of this solution.Attempts to verify the KZ solution heavily focus on the scaling exponent γ [7][8][9][10], with only a handful of them targeting the Kolmogorov constant C.Among the latter, recent experimental validations [e.g., 11,12] lead to larger discrepancies with theory (as much as a factor of 20 difference), partly due to the challenge of precisely measuring P from experimental data.The studies that are relatively successful consider well-controlled numerical simulations, including turbulence of capillary waves [13][14][15] and Bose-Einstein Condensates [16], which provides values of C respectively about 40% above and 10% below the corresponding theoretical values.These two results signify that for a given power-law spectrum with theoretical slope γ, the energy flux computed from simulations are respectively about 0.5 and 1.3 (≈ (1/1.40) 2 , 1.1 3 ) times of that from the KZ solution.* ahrabski@umich.edu The scarce and limited success in verifying the Kolmogorov constant indicates an insufficient understanding of the validity of the WKE for stationary WT.This issue is more subtle than verification of WKE at an evolving state [e.g., 17] where spectral evolution serves as a natural measure of success.We are mainly interested in two fundamental questions in this work: (1) What are the major obstacles in obtaining theoretical value of C in simulations of dynamical equations?In particular, what does it take to bring the numerical value of C further closer to the theoretical value than in previous validations?(2) How is KZ spectrum (and more generally, the WT closure) realized in one-dimensional (1D) systems?We note that 1D WT is in a sense more difficult to describe than in higher dimensions [16][17][18] due to far fewer interactions for any given spectral range.One example of a 1D system is the Majda-McLaughlin-Tabak model [7], for which the theoretical value of γ has been notoriously difficult to reproduce numerically for more than 20 years.This puzzle was resolved only recently in [9], where it was shown that the theoretical value of γ can be recovered with a much wider inertial range than those in previous studies.The Kolmogorov constant for this 1D system, on the other hand, has never been numerically studied.
With the two above questions in mind, we perform a numerical study of the MMT model focusing on the Kolmogorov constant in the stationary state.One feature distinguishing our current study from all previous studies is that we numerically probe the kinetic limit, by weakening nonlinearity while making the domain large.We show that the theoretical value of C can be recovered as a result of two limiting processes: (1) as we approach the kinetic limit, the dynamics of the MMT model converge to the WKE description, evaluated through the closure model and energy flux; (2) as we enlarge the inertial range, the value of C computed from the collision integral converges to the theoretical value to very high precision.In component (1), we find that quasi-resonances can lead to difference between MMT simulation and WKE prediction when the former is taken outside the kinetic limit as in [16].In component (2), we find that an inertial range of at least 3.5 decades is needed for the convergence of C.
The MMT model is a family of nonlinear dispersive equations for a complex field ψ(x) = ψ ∈ C that are widely used in the study of WT [7,[19][20][21].The MMT equation of interest to this work reads where the derivative operator |∂ x | 1/2 produces a dispersion relation ω k = |k| 1/2 , the same as surface gravity waves.We consider the MMT equation on a 1D periodic domain of length L, where we have ψ(t) = k∈ΛL ψk (t)e ikx , with Λ L ≡ 2πZ/L.The statistical description of (1) begins by defining the wave action spectrum , where the • denotes an ensemble average (or a time average for statistically stationary data).In deriving the WKE governing n k , a statistical closure needs to be taken which connects 4th-order correlation of ψk to n k , in the form of where Ω denotes the frequency mismatch of the four wave modes k 1 , k 2 , k 3 and k.The closure (2) can be derived in various ways, assuming quasi-Gaussian statistics [1], or more recently under the less restrictive assumption of a field with random phases and amplitudes [22][23][24].
Depending on different methods of derivation, f (Ω) takes different forms of a broadened delta function, namely a Lorentzian form [25] or sinc-like functions [26,27].At the kinetic limit, we have f (Ω) → δ(Ω) and reach the WKE.
One can further seek stationary solutions to (3) using the so-called Zakharov transformation [25].Of interest here is the KZ solution associated with a finite energy flux from small to large k, taking the form of n k = CP 1/3 k −1 .Included in the supplemental material [28] is a full derivation of the KZ spectrum, including a correction to the previous work [19] leading to a new value of C = 0.2984.

III. METHODS OF NUMERICAL STUDY
We simulate (1) via the pseudospectral method developed in [7].To create a steady forward energy cascade, we add Gaussian forcing in the range of 10 ≤ |k| ≤ 20 and a dissipation term −iν k ψk to the RHS of (1), with ν given by For each simulation, we start from a quiescent field and simulate until a stationary spectrum is reached, with the final nonlinearity level measured by ǫ = H 4 /H 2 , where H 2 and H 4 are the linear and nonlinear components of total energy [7], respectively.To numerically approximate the kinetic limit, we choose four forcing strengths leading to four final nonlinearity levels, and for each forcing strength, we conduct simulations on domain of sizes L ∈ [2π, 4π, 8π, 16π].As L increases from 2π to 16π, the resolution in k space progressively doubles with ∆k decreasing from 1 to 1/8, with k max = 1024 kept for all simulations.In doing so extra care needs to be taken to scale the forcing and dissipation parameters with the domain size L, in order to keep quantities such as ǫ and total energy nearly constant across different domain sizes.Details on how this is done are included in the supplemental materials [28].
With numerical data available, we are interested in studying the behavior of the closure model (2), especially in the context of a large number of quartets forming the energy flux.The basis for this analysis is an exact evaluation and decomposition of energy flux developed in [29]: where k b is the wave number through which the time- The measured closure function f (Ω) (--) for each tested L = 16π case, denoted from highest to lowest ε by green, magenta, red, and blue, respectively.A fitted Lorentzian closure f (Ω) for each case (----).
Eq. ( 5) describes the decomposition of energy flux into contributions from quartets with the wavenumber condition satisfied and with a mismatch of Ω in frequency condition.
The formulation of P Ω in ( 6) can be derived directly from (1) without any assumption [29].With the closure model ( 2) substituted in ( 6), we have By computing the LHS of ( 7) through ( 6) and the RHS via n k taken from our simulation data, we can directly compute the functional form of f (Ω).Comparison of the numerically-resolved f (Ω) with the analytical function then provides us a metric to evaluate the validity of the closure model.We note that this method evaluates the performance of (2) over a large number of quartets, which is shown in [29] as the only meaningful way to study the closure model with a single simulation.In addition, since the RHS of ( 7) is exactly the energy flux calculated from a discrete form of WKE with a broadened delta function, named quasi-resonant WKE (QRWKE) in [14], the closeness between numerical and analytical f (Ω) also indicates the accuracy of the WKE (or QRWKE) in reproducing the energy flux in dynamical simulations.

IV. RESULTS
We begin by showing in figure 1 the stationary spectra for all 16 simulations varying nonlinearity level ǫ ∈ (0.0066, 0.067) and domain size L. At high nonlinearity levels, we see good agreement in spectral form for all L, indicating that even the smallest domain size L = 2π is sufficient to capture the large-L dynamics.At low nonlinearity (especially ǫ = 0.0066), however, we see that the spectrum varies substantially as L changes.In particular, the secondary peaks occurring at small L reflect finitesize effects.They disappear as L increases, suggesting a transition from the discrete to the kinetic turbulence regime [8,30].
Also shown in figure 1 are the KZ solutions, with P computed though k b = 300 for L = 16π for each ε.While we see the spectral slope gets closer to the KZ value of γ = −1 as as ε decreases, the inertial interval also shrinks, and perhaps even departs slightly from a true power-law.Additionally, the spectral level of the numerical solution does not get closer to the KZ solution as nonlinearity is decreased, indicating that the numerically-resolved Kolmogorov constant does not agree better with its theoretical value.We note that this is not in contradiction with the major theme of the paper.As discussed below, we need to consider two limiting processes to precisely reproduce the theoretical value of C: one taking the kinetic limit, and the other increasing the width of the inertial range.
Before discussing the two limiting processes in detail, we would like to briefly mention another set of important results that are made available due to our detailed analysis.As demonstrated in the supplemental material [28], the energy flux follows a Gaussian distribution, whereas the energy dissipation rate follows a log-normal distribution.Both have a standard deviation that scales with 1/ √ L. The log-normal distribution on the dissipation rate, first observed in wave turbulence here, is especially interesting as it has also been described in flow turbulence [31,32] with a sound theoretical explanation yet to be developed.

A. Limiting Process 1: Kinetic Limit
In this section, we will show that as the kinetic limit is taken, the measured energy flux indeed converges to the prediction of the WKE.Following methods introduced in §III, we plot f (Ω) for different nonlinearity levels at L = 16π, together with the analytical Lorentzian form a/π(a 2 + Ω 2 ) [25] with values of a that best fit the data.We see that the analytical form fits the data remarkably well, and that as nonlinearity is decreased with sufficient domain size, the function f (Ω) approaches a true delta function.These results indicate that the long-time dynamics indeed follows the WT closure and converges to WKE description as the kinetic limit is taken.Due to the broadening of f (Ω), these results also suggest that dynamics away from the kinetic limit may be better represented by the QRWKE.Since the Lorentzian form of the delta functions contains long tails, it is expected that a large number of quasi-resonances are active in dynamic simulations away from the kinetic limit.These quasiresonances are therefore key factors resulting in difference in energy flux from dynamical simulations and WKE calculations as seen in [16].

B. Limiting process 2: Wide Inertial Range
The convergence of dynamics to the WKE prediction at the kinetic limit does not guarantee that the simulated spectrum converges to the KZ solution, as we have demonstrated.As a result, the Kolmogorov constants evaluated for the simulated spectrum using the dynamic flux (6) and kinetic flux (7) with f (Ω) → δ(Ω) are very close to each other but away from the theoretical value (see symbols in figure 3b).To understand this situation, recall that, given a valid WKE, the realization of the KZ solution requires the dominance of local interactions, which in turn requires a sufficiently wide inertial range.The width of inertial range seems to be especially important for 1D models such as MMT, as evidenced in [9] on the sensitivity of γ to the width.It is our objective to show next that our numerical solution is indeed contaminated by non-local interactions, and that a much wider inertial range is needed to precisely recover the theoretical value of C.
We first quantify the influence of forcing peak (as the source of non-local interactions) to the high-wavenumber portion of spectrum.Considering the stationary state, we explicitly decompose the average small-scale dissipation rate via P d = J nl + J l , where J is the small-scale energy evolution due to nonlinear interactions with subscripts nl and l denoting contributions from non-local and local interactions.An interaction is considered nonlocal if at least one wave number resides in forcing range (10 ≤ |k| ≤ 20).Refining methods developed in [33] (see details in [28]), we explicitly compute J nl and plot the ratio J nl /P d in figure 3a.At low nonlinearities, about 40% of the inter-scale energy flux is due to non-local interactions, clearly violating the assumptions of the KZ spectrum.Interestingly, the non-local contribution decreases as nonlinearity is increased, explaining the longer power-law range at higher nonlinearity seen in figure 1.
We finally demonstrate precise convergence to the theoretical value of C as the width of inertial range is made larger.This is achieved by evaluating kinetic flux P W KE for an idealized spectrum n k = k −1 with cutoffs at both low and high wavenumbers k a and k c .Here ka ω k ∂n k /∂t with ∂n k /∂t computed from the collision integral of the WKE (3), which is exactly RHS of (7) with f (Ω) → δ(Ω).Specifically, we keep k a = 10 and progressively move k c to higher wavenumber to represent the widening of the inertial range, with P W KE evaluated at k b = 560.Figure 3b plots the value of C as a function of k c /k a .At k c /k a ≈ 100, some discrepancy is seen between C evaluated on the simulated and idealized spectra, due the differences in form between the two spectra.With the increase of k c /k a , we see that value of C converges to the theoretical value with very high precision.We find that an inertial range of about 3.5 decade is needed to obtain converged and accurate result for C.This is much longer than what is needed to recover the theoretical γ in MMT model [9] and what is needed for C in high-dimensional models [16].We note that for the MMT model we consider, γ = −1 is far from the boundaries of the locality window −7/4 < γ < 1/2 [19].The long inertial range needed to recover the KZ solution is therefore more likely a general feature of the 1D wave turbulence.

V. CONCLUSION
We present a detailed numerical study on the MMT model, focusing on the realization of the Kolmogorov con-stant and the WT closure.We show that the Kolmogorov constant can be precisely recovered following two limiting processes: the first as kinetic limit is taken where the dynamics converges to the WKE prediction in terms of the WT closure model realized through energy flux; the second as the inertial-range spectrum is made sufficiently wide so that enough interactions over the inertial range are captured.
The derivation of the KZ spectrum begins with the Wave Kinetic Equation (WKE), provided again for convenience: While it isn't in general necessary [1], we adopt the popular assumption of an isotropic spectrum.In one-dimensional systems, this takes the form of n k = n −k .Each of k 1 , k 2 , and k 3 can take a positive or negative sign in the delta function of (1), producing 8 possibilities.However, not every combination of these signs produces a non-trivial resonance.A trivial resonance is one that satisfies the resonance conditions by having k 1 = k 3 and k 2 = k, or k 2 = k 3 and k 1 = k.For these cases, the integrand of (1) takes a zero value, so they may be ignored.After these trivial cases are removed, we are left with It will later become important to integrate over the resonant manifold (described by the δ-functions), which is much easier to interpret as quadratic functions in ω = k 1/2 .Therefore, we next rewrite the above equation as an integral over ω, and replace (on the LHS) n k , the spectral density in k, with N ω , the spectral density in ω.This leads to where n ω = n(k(ω)), and we have used the facts that dk = 2ωdω and N (ω)dω = n(k)dk + n(−k)dk = 2n(k)dk.We note that the factor of 2 on the spectral element relation is necessary for a consistent and correct flux definition.We now assume that n ω = Aω γ , where γ refers to an arbitrary exponent that will later be used to determine the KZ exponents.Substituting this into (3), we are left with Next we employ the Zakharov transformations [1,2], which are a set of conformal transformations one applies to the integrand that result in the reduction of the sum of delta functions to a single delta function.This new structure of the integrand will allow us to (a) explicitly see the zeros of the equation and (b) explicitly compute the Kolmogorov constant C. See [3] for an intuitive, geometric description of how these transformations achieve these effects assuming only a self-similar spectrum.We distribute the sum of delta functions to expand the integrand into 4 terms, and we handle each separately.The first term, is the form onto which we will map the other terms.We have used a superscript (1) to denote we are referring to the first term.Now, as an example, we manipulate the second term in full, to which we apply the following transformations: , and ( 6) becomes, after some reduction, We note that the identity δ(ax)dx = δ(x)/|a|dx is used to simply the above expression.The integrand of ( 7) almost reflects ( 5) with an additional factor of (ω/ω ′ 1 ) 3γ+5 .If one carefully looks at the signed terms in the equation, however, it becomes apparent that certain indices have become switched as a result of our transformation.The choice of indices is arbitrary, so we renumber them according to 123 → 132.This leaves us with In this form, the symmetry with (5) is obvious.We perform the remaining Zakharov transformations (see [2] for the forms of the other transformations), sum the four terms, and drop the primes from our notation.This results in with y ≡ −3γ −5.Careful consideration of the resonance conditions reveals that ( 9) is an integral over the intersection of the plane ω 1 + ω 2 − ω 3 − ω = 0 and the sphere ω 2 1 + ω 2 2 + ω 2 3 = ω 2 for any given ω.Thus, we do not need to consider integrating over any ω i > ω.This enables a reparameterization in terms of some This is the form of the WKE that allows for the solution of γ corresponding to stationary solutions.One can see by the second product in the integrand that if γ = 0, the RHS of ( 10) is identically 0. This corresponds to equipartition of wave action, i.e., n ω is constant.Also, if γ = −1, then the second product is identically 0 whenever the resonance condition is satisfied, also leading to a 0 of the collision integral.This solution n ω = A/ω corresponds to the Rayleigh-Jeans spectrum.Neither of these solutions correspond to wave turbulence, but rather are equilibrium solutions.The KZ solutions are given by y = 0 and y = 1, which produce 0's of the collision integral by the same arguments as the equilibrium solutions.Setting y = 0, one obtains n ω = Aω −5/3 , and for y = 1, one obtains n ω = Aω −2 .To determine which of these out-of-equilibrium spectra correspond to the forward cascade of energy and which corresponds to the inverse cascade of wave action, one may use the ordering of the exponents γ relative to the equilibrium spectra [1], or Zakharov's method via evaluating (10) [4], while in both cases being careful to ensure the flux directions are not non-physical via comparison to the equilibrium spectra [1,2,4].In [4], it is demonstrated that the forward and inverse cascade for our system have physical cascade directions and that y = 1 (with γ = −2) corresponds to the forward cascade.
Next, we derive the Kolmogorov constant C. The energy flux through frequency ω is defined by a control volume argument in spectral space as where an negative sign is introduced to ensure that a positive flux corresponds to a cascade of energy from large to small scales (long to short time scales via the dispersion relation).For the forward cascade with y = 1, the computation of P (ω) involves the limit of an indeterminate quantity, which can be obtained via L'Hospital's rule to be The limit of the desired derivative, S, is given by where we have used the fact that lim y→1 dx y dy = x ln x.In the next section, we develop a precise method to numerically evaluate (13) to find S = 0.09353.We now compute the relationship between A and P (ω), revealing both the Kolmogorov constant C via resulting in C = 0.2984.Given that C is positive, we can now be sure that the cascade direction is correct.Thus, the KZ spectrum associated with the forward cascade process in our MMT equation is given by n ω = 0.2984P 1/3 ω −2 , or, via the linear dispersion relation, n k = 0.2984P 1/3 k −1 .We note that this value of C is different from the result of Zakharov et al. [4], due to their missing factor of 2 in the spectral element relation N (ω)dω = 2n(k)dk.

II. NUMERICAL INTEGRATION OVER THE RESONANT MANIFOLD
A. Evaluation of S We are interested in integrating (13).For simplicity, we will refer the non-delta part of the integrand by f (ξ 1 , ξ 2 , ξ 3 ) so that This leaves The first of these delta functions is linear in ξ 1 .Making use of the property we can integrate over ξ 3 to obtain Integrating over the region ∆(ξ 1 , ξ 2 ) ≡ {0 < ξ 1 < 1, 0 < ξ 2 < 1, 1 < ξ 1 + ξ 2 < 2} simply ensures 0 < ξ 3 < 1.We would like to now apply (17) again, however we require a transformation so that the argument ξ 2 1 +ξ 2 2 +ξ 2 3 −1 is of the required form.To do this, we transform the inner integral to one with respect to du, where u Now, we may apply (17), being careful to include only the part of u(ξ 1 , ξ 2 ) = 0 that lies in ∆(ξ 1 , ξ 2 ).After some manipulation of our definition of u we find that, of the two branches for which u = 0, the one with is in the region ∆(ξ 1 , ξ 2 ).After applying (17), This form of S is suitable for numerical integration, where the integrand as ξ 2 → 1 (from below) can be shown to approach zero via L'Hospital's rule.

B. Evaluation of the WKE Collision Integral
Just as in the evaluation of dI(y) dy , this explicit algebraic approach works to evaluate (3).To start, we rewrite (3) in terms of the function f For each δ-function with quadratic arguments in ω, we perform the procedure outlined in the previous subsection.After a good deal of simplification, this results in the following expression for ∂Nω ∂t : (u),ω (u),ω 2 (u),ω (u),ω (u),ω The first delta function on the RHS of ( 22) corresponds to ∂N ω /∂t, the next two delta functions (due to symmetry) produce identical contributions and are combine in ∂N This is a form of the collision integral that is suitable for numerical integration.When computing P W KE directly from the collision integral in the main text, n k ≡ n(k(ω)) is evaluated explicitly from the idealized, truncated KZ spectrum given by When the simulated spectrum is instead used in this evaluation (for the two points in figure 3b of the main text), we use piecewise-linear interpolation of the simulated spectrum to evaluate n k for k / ∈ Λ L .In both of these cases, n k = 0 after some large value of k with a corresponding u, meaning that the upper integral bounds of ∂N (3a) ω /∂t and ∂N (3b) ω /∂t in our numerical evaluations can be taken to be finite.All integrals are evaluated via adaptive quadrature.

III. SELECTION OF NUMERICAL PARAMETERS
For a fair comparison between domains of different L, we would like to keep key quantities like the length-averaged Hamiltonian density H, ε, and P approximately constant as L increases, while allowing other quantities to vary.Choosing the parameters of each numerical simulation to achieve this goal turns out to be a non-trivial task.Given the definition of our Fourier series the spectral amplitude ψk is normalized such that the length-averaged total action does not depend on L, as can be seen via Parseval's theorem [1], The key parameters that control the energy content of our simulations are forcing and dissipation.Thus, we attempt to keep the action injection rate, controlled by the standard deviation of the Gaussian forcing σ F , and the dissipation

FIG. 3 .
FIG. 3. (a) The relative contribution of (non-local) interactions with the forcing range to the average small-scale dissipation rate for varying ε, plotted for the (representative) L = 16π case.(b) Value of C for varying inertial interval length kc/ka.( ) denotes C computed via the kinetic flux for the simulated spectrum.(•) denotes C computed via dynamic flux for the simulated spectrum.