1. Introduction
Wave turbulence theory describes the statistical properties of ensembles of weakly nonlinear interacting waves, with rich applications in many physical contexts, e.g. ocean waves (Zakharov & Filonenko Reference Zakharov and Filonenko1967; Nazarenko & Lukaschuk Reference Nazarenko and Lukaschuk2016), acoustics (L'vov et al. Reference L'vov, L'vov, Newell and Zakharov1997), magnetohydrodynamics (Galtier et al. Reference Galtier, Nazarenko, Newell and Pouquet2000), quantum turbulence (Nazarenko & Onorato Reference Nazarenko and Onorato2006) and others. The centrepiece of wave turbulence theory is a wave kinetic equation (WKE), which describes the time evolution of the wave action spectrum as an integral over wave–wave interactions. The WKE yields a stationary Kolmogorov–Zakharov power-law solution associated with a constant flux, which has been observed in many wave systems.
Quantitative validations of the WKE and its predictions in different physical contexts have been a prominent topic in the wave turbulence community for decades. These studies include extensive numerical and experimental validations of the spectral slope and energy flux (or Kolmogorov constant) of the Kolmogorov–Zakharov solutions (e.g. Pan & Yue Reference Pan and Yue2014; Falcon & Mordant Reference Falcon and Mordant2022; Hrabski & Pan Reference Hrabski and Pan2022; Zhang & Pan Reference Zhang and Pan2022a; Zhu et al. Reference Zhu, Semisalov, Krstulovic and Nazarenko2023), numerical validation of the initial spectral evolution predicted by the WKE (e.g. Banks et al. Reference Banks, Buckmaster, Korotkevich, Kovači and Shatah2022; Zhu et al. Reference Zhu, Semisalov, Krstulovic and Nazarenko2022) and rigorous mathematical justification of the WKE (e.g. Buckmaster et al. Reference Buckmaster, Germain, Hani and Shatah2021; Deng & Hani Reference Deng and Hani2021, Reference Deng and Hani2023). Although successes in verifying the WKE have been reported in many cases, situations where WKE predictions fail have also been identified, such as circumstances associated with finite-size effects (e.g. L'vov & Nazarenko Reference L'vov and Nazarenko2010; Hrabski & Pan Reference Hrabski and Pan2020; Zhang & Pan Reference Zhang and Pan2022b), coherent structures in the field (e.g. Rumpf, Newell & Zakharov Reference Rumpf, Newell and Zakharov2009) and the strong turbulence regime (e.g. Chibbaro, De Lillo & Onoroto Reference Chibbaro, De Lillo and Onoroto2017).
In spite of the significant advancement in understanding the WKE, there exists a series of relevant studies on surface gravity waves that remain largely unexplained. A representative work of this series is the paper by Annenkov & Shrira (Reference Annenkov and Shrira2009), which considers the spectral evolution when a stationary gravity wave field is subjected to an abrupt perturbation of external wind forcing. According to the WKE of gravity waves, in the form of 
$\partial n/\partial t \sim n^3$ with 
$n\sim \varepsilon ^2$ the wave action spectrum, one would expect an evolution with the kinetic time scale of 
$O(\varepsilon ^{-4})$. However, simulations of the dynamic equation, specifically the Zakharov equation (Zakharov Reference Zakharov1968), show a faster evolution with a dynamic time scale of 
$O(\varepsilon ^{-2})$. We remark that this ‘fast’ spectral evolution refers to the scaling of time scale with 
$\varepsilon$, instead of the evolution rate in its absolute value. Although with some uncertainties in measurement, the ‘fast’ spectral evolution is also observed in wave-tank experiments (e.g. Autard Reference Autard1995; Waseda, Toba & Tulin Reference Waseda, Toba and Tulin2001) and field studies (van Vledder & Holthuijsen Reference van Vledder and Holthuijsen1993) when the wind exhibits a sudden change in speed or direction, as well as numerical simulations for the initial evolution of some wave spectra (Dyachenko et al. Reference Dyachenko, Trulsen, Krogstad and Socquet-Juglard2003). In a study by Annenkov & Shrira (Reference Annenkov and Shrira2018), an attempt is made in explaining the fast spectral evolution using the so-called generalised WKE, but only limited success is achieved. This unexplained issue is concerning since the WKE of surface gravity waves (also known as Hasselmann's kinetic equation; Hasselmann Reference Hasselmann1962) is currently used in modern wave modelling codes whose reliability is pertinent to weather forecasting, climate modelling and navigation (Janssen Reference Janssen2004).
In this paper, we revisit the spectral evolution problem in the context of the one-dimensional (1-D) Majda–McLaughlin–Tabak (MMT) equation with gravity wave dispersion relation. The MMT model is favourable in the sense that it captures the essential dynamics of wave turbulence while being exempt from the complexities associated with surface gravity waves. One of these critical complexities is the existence of quasi-resonant three-wave interactions in surface gravity waves (an issue under investigation in the field), which can affect the observed spectral evolution time scale. The mathematical procedure to remove the quadratic terms in deriving the WKE of surface gravity waves is also not fully justified for finite wave steepness 
$\varepsilon$. In addition, with a 1-D model, we are able to perform large numbers of ensemble simulations, with which the statistical properties can be reliably computed through ensemble averages. We conduct the study and analysis closely following Annenkov & Shrira (Reference Annenkov and Shrira2009), i.e. with numerical set-ups and evaluation of properties as consistent as possible, except using the 1-D MMT model.
For a given stationary wave field, we show that evolution with kinetic and dynamic time scales can be observed for weaker and stronger forcing perturbations, respectively. The transition from the former to the latter regime corresponds to a forcing perturbation that triggers a spectral evolution time scale comparable to the linear time scale, i.e. the wave period of some representative mode, violating the basis of the WKE. We further show, through a study varying the energy of the base stationary wave field, that this violation of time-scale separation is primarily a result of the spectrum far from the stationary state after a sufficiently strong forcing perturbation is applied, rather than the increase of overall nonlinearity level of the wave field. We finally develop a set-based interaction analysis, which allows us to understand that the spectral evolution for modes in the inertial range is predominantly governed by their (mostly non-local) interactions with the low-wavenumber condensate (or regions sufficiently filled with energy) induced by the forcing perturbation.
2. Methodology
2.1. The MMT model
We study the evolution of random wave fields through the 1-D MMT equation (Majda, McLaughlin & Tabak Reference Majda, McLaughlin and Tabak1997), which is a family of nonlinear dispersive wave equations widely used to study wave turbulence problems (e.g. Cai et al. Reference Cai, Majda, McLaughlin and Tabak1999; Zakharov, Pushkarev & Dias Reference Zakharov, Pushkarev and Dias2001; Chibbaro et al. Reference Chibbaro, De Lillo and Onoroto2017; Hrabski & Pan Reference Hrabski and Pan2022):
\begin{equation} \textrm{i}\frac{\partial\psi}{\partial t}={\lvert{\partial _x}\rvert}^\alpha \psi+ \lambda{\lvert{\partial _x}\rvert}^{\beta/4}({\lvert{\lvert{\partial _x}\rvert}^{\beta/4}\psi\rvert}^2 {\lvert{\partial _x}\rvert}^{\beta/4}\psi), \end{equation}
where 
$\psi (x,t)$ is a field taking complex values. The parameter 
$\beta$ controls the nonlinearity formulation and 
$\alpha$ controls the dispersion relation 
$\omega (k)=|k|^\alpha$ with 
$\omega$ the frequency and 
$k$ the wavenumber. Here 
$\lambda =-1$ and 
$1$ represent the focusing (defocusing) nonlinearity, respectively (Cai et al. Reference Cai, Majda, McLaughlin and Tabak2001). The MMT equation (2.1) conserves both the total Hamiltonian and wave action. In our study, we use 
$\alpha =1/2$ to mimic the dispersion of surface gravity waves, and use 
$\beta =0$ for convenience which is consistent with a portion of the original study of the MMT equation (Majda et al. Reference Majda, McLaughlin and Tabak1997) and is widely used in other studies (e.g. Cai et al. Reference Cai, Majda, McLaughlin and Tabak1999; Rumpf & Newell Reference Rumpf and Newell2013; Rumpf & Sheffield Reference Rumpf and Sheffield2015). We note that although 
$\beta =0$ is not entirely representative of surface gravity waves, the primary conclusion of the paper does not depend on a specific choice of 
$\beta$. We present the results for the defocusing case (
$\lambda =1$) in the main paper, and those for the focusing case in Appendix B, with both cases exhibiting similar physics regarding the spectral evolution time scales.
From a standard wave turbulence consideration, a statistical description of the wave field can be obtained by defining the wave action spectrum 
$n_{k}=\langle {\lvert {\hat {\psi }_k}\rvert }^2\rangle$, with 
$\hat {\psi }_k$ the Fourier transform of 
$\psi$ and the angle brackets denoting an ensemble average. Under conditions of weak nonlinearity, random phases and infinite domain, the time evolution of 
$n_k$ is governed by the WKE (for 
$\beta =0$):
\begin{align} \frac{\partial n_{k}}{\partial t} &= 4{\rm \pi}\int (n_{1}n_{2}n_{3}+n_{1}n_{2} n_{k}-n_{1}n_{3}n_{k}-n_{2}n_{3}n_{k}) \nonumber\\ &\quad \times\delta(\omega_{1}+\omega_{2}-\omega_{3}-\omega) \delta(k_{1}+k_{2}-k_{3}-k)\,\textrm{d}k_{1}\,\textrm{d}k_{2}\,\textrm{d}k_{3}, \end{align}
where 
$\delta$ is the Dirac delta function. According to (2.2), the spectrum evolves with a kinetic time scale, i.e. 
$n_k$ experiences significant change with a time scale of 
$O(\varepsilon ^{-4})$, and 
$\partial n_k/\partial t\sim O(\varepsilon ^6)$ over the evolution, with 
$\varepsilon \sim \sqrt {n_k}$ a measure of wave steepness. Hereafter, we refer to these relations as kinetic scaling in this paper, which is in contrast to the dynamic scaling of a time scale of 
$O(\varepsilon ^{-2})$ and 
$\partial n_k/\partial t \sim O(\varepsilon ^4)$ that one can directly obtain from (2.1).
2.2. Numerical procedure
We simulate (2.1) with 4096 modes (before de-aliasing) on a periodic domain of size 
$L=2{\rm \pi}$, with the addition of forcing and dissipation terms. The forcing is in white-noise form, given by
\begin{equation} F =\begin{cases} F_{r}+\textrm{i}F_{i}, & 4\le k \le 13,\\ 0, & \text{otherwise}, \end{cases}\end{equation}
with 
$F_r$ and 
$F_i$ independently drawn from a Gaussian distribution 
$\mathcal {N}(0, \sigma ^2)$. The dissipation is introduced with the addition of two hyperviscosity terms
\begin{equation} \left.\begin{aligned} D_{1} =\begin{cases} -\textrm{i}\nu_{1}\hat{\psi}_{k}, & k\ge 900,\\ 0, & \text{otherwise}, \end{cases} \\ D_{2} =\begin{cases} -\textrm{i}\nu_{2}\hat{\psi}_{k}, & k\le 4,\\ 0, & \text{otherwise}, \end{cases}\end{aligned}\right\} \end{equation}
at small and large scales, respectively. Since the MMT model supports the inverse cascade, it is necessary to use large-scale dissipation to avoid energy accumulation at large scales. The dissipation coefficients are fixed to be 
$\nu _{1}=1\times 10^{-14}(k-900)^8$ and 
$\nu _{2}=3k^{-4}$ for all numerical experiments. The numerical schemes used for the simulation are discussed in detail in previous papers (Hrabski & Pan Reference Hrabski and Pan2020, Reference Hrabski and Pan2022).
In order to reproduce the physical scenario as in Annenkov & Shrira (Reference Annenkov and Shrira2009), we perform the simulation with two stages. In stage 1, we simulate (2.1) with (2.3) and (2.4) using 
$\sigma =\sigma _0$ for sufficient time to reach a stationary wave field (and spectrum), starting from an initial spectrum that exponentially decays in 
$|k|$. In stage 2, we add a perturbation 
${\rm \Delta} \sigma$ to the forcing magnitude (leading to 
$\sigma =\sigma _0+{\rm \Delta} \sigma$) and study how the spectrum evolves subject to the perturbation. In both stages, we use a time step 
${\rm \Delta} t=0.012$. For the main cases studied in this paper, we consider the stationary state obtained with forcing 
$\sigma _0=0.004$ in stage 1, and consider forcing perturbations 
${\rm \Delta} \sigma \in [0.033,1.461]$ in stage 2 that is sufficient to cover the physical regimes of our interest. The process is simulated with an ensemble of 65 000 simulations with different random seeds in the forcing, and all statistical results related to 
$n_k$ presented below are obtained directly from the average over the whole or part of the ensemble that is statistically convergent.
3. Results
Figure 1 shows the evolution of a spectrum in a typical case with 
$\sigma _0=0.004$ and 
${\rm \Delta} \sigma = 0.033$. We see that a stationary power-law spectrum forms at the end of stage 1, which serves as the base state of the problem. As a forcing perturbation is excited in stage 2, the forcing scales (as well as scales slightly larger than that) first experience an abrupt growth, forming a condensate at large scales. The growth is then propagated to smaller scales, which eventually fill in the full spectral range and form the final stationary power-law spectrum. We note that such behaviour regarding the propagation of growth is closely related to the fact that (2.1) under selected values of 
$\alpha$ and 
$\beta$ is an infinite-capacity system (Nazarenko Reference Nazarenko2011; Newell & Rumpf Reference Newell and Rumpf2011). Hereafter, for simplicity, we assign time 
$t=0$ to the stationary state before a forcing perturbation is applied, and 
$t=t_f$ to the time when the final stationary state is formed. We measure the wave steepness of each state as 
$\varepsilon =\sqrt {E_L}$ with 
$E_L$ the linear energy (or Hamiltonian) of the system. Accordingly, we use 
$\varepsilon _0$ and 
$\varepsilon _f$ as the wave steepness at 
$t=0$ and 
$t=t_f$. We note that this definition of 
$\varepsilon$ is consistent with that in Annenkov & Shrira (Reference Annenkov and Shrira2009) with the additional factor of constant peak wavenumber 
$k_p$ omitted in our definition.
Figure 1. A typical case of the evolution of wave action spectrum, where 
$n(k,0)$ (solid blue) and 
$n(k,t_f)$ (solid red) represent the stationary spectra before the forcing perturbation (stage 1) and after sufficient evolution with forcing perturbation applied (stage 2). The intermediate spectra in the evolution between 
$n(k,0)$ and 
$n(k,t_f)$ are plotted (dashed red).
We are interested in the growth rate in the transient period in stage 2, based on which we can evaluate 
$\partial n_k/\partial t$ and the associated scaling with wave steepness. For this purpose, we follow Annenkov & Shrira (Reference Annenkov and Shrira2009) to measure 
$\partial n_k/\partial t$ for a given mode 
$k$ using data over the interval of 5 %–40 % of the full modal growth. More specifically, if we define a normalised wave action spectrum 
$\tilde {n}(k,t)=(n(k,t)-n(k,0))/(n(k,t_f)-n(k,0))$, we can then measure 
$\partial \tilde {n}_k/\partial t$ over the range of 
$\tilde {n}_k \in [5\,\%, 40\,\%]$ via a least-squares fit, and then scale back to compute 
$\partial n_k/\partial t$. Figure 2 shows the evolution of 
$\tilde {n}(k,t)$ for three selected modes in the inertial range, as well as the growth rate evaluated using the 
$[5\,\%, 40\,\%]$ interval. We note that each modal growth exhibits a short exponential regime followed by an algebraic regime, and that our 
$[5\,\%, 40\,\%]$ covers the growth predominantly in the algebraic regime. Furthermore, we have tested that the obtained results are not sensitive to a variation of bounds of the interval to [25 %, 60 %] and likely more. Therefore, the evaluated growth rate over the chosen interval 
$[5\,\%, 40\,\%]$ is robust to capture the typical spectral growth after the forcing perturbation is imposed.
Figure 2. A typical case of the evolution of modal wave action at three wavenumbers 
$k=150$ (solid green), 
$k=500$ (solid blue) and 
$k=800$ (solid red) in the inertial range. The growth rates 
$\partial \tilde {n}_k/\partial t$ measured over the range of 
$\tilde {n}_k \in [5\,\%, 40\,\%]$ are indicated (dashed).
3.1. Scaling of spectral growth rate with 
$\varepsilon$
To evaluate the scaling of 
$\partial n_k/\partial t$ with the wave steepness 
$\varepsilon$, we use the final-state 
$\varepsilon =\varepsilon _f$ as a representative value of the wave steepness for the case (cf. Annenkov & Shrira Reference Annenkov and Shrira2009). Figure 3 shows 
$\partial n_k/\partial t$ of three modes in the inertial range for a broad range of 
$\varepsilon \in [0.022,0.124]$, obtained from 22 cases with 
${\rm \Delta} \sigma \in [0.033,1.461]$ starting from a base state with 
$\varepsilon _0=0.010$ resulting from 
$\sigma _0=0.004$. As we see in figure 3, the kinetic scaling 
$\partial n_k/\partial t \sim O(\varepsilon ^{6})$ is realised in the range of 
$\varepsilon \lesssim 0.05$, corresponding to small 
${\rm \Delta} \sigma$ with 
${\rm \Delta} \sigma \lesssim 0.36$. For larger forcing perturbations, we find dynamic scaling in the range of 
$\varepsilon \gtrsim 0.05$. We remark that in the previous work for gravity waves (Annenkov & Shrira Reference Annenkov and Shrira2009), only dynamic scaling is observed. This is possibly due to those studies only being conducted for strong forcing perturbations, or other reasons that we will leave for future study.
Figure 3. Rate of spectral evolution 
$\partial n_k/\partial t$ as a function of nonlinearity 
$\varepsilon$ for (a) 
$k=150$ (green triangles), (b) 
$k=500$ (blue triangles) and (c) 
$k=800$ (red triangles), with the kinetic 
$\varepsilon ^{6}$ and dynamic 
$\varepsilon ^{4}$ scalings indicated by the dashed lines.
We next investigate why the ‘fast’ dynamic scaling becomes relevant for the spectral evolution, a critical question that is not answered in previous works (e.g. Annenkov & Shrira Reference Annenkov and Shrira2006, Reference Annenkov and Shrira2009, Reference Annenkov and Shrira2018). The realisation of dynamic scaling indicates that the WKE (2.2) (or more generally, wave turbulence theory) must break down. One situation in which this could happen is a violation of time-scale separation, i.e. when the nonlinear modal time scale becomes comparable to the linear modal time scale. In particular, the linear time scale of a mode is given by the modal wave period:
\begin{equation} \tau_{lin}=\frac{2{\rm \pi}}{\omega}. \end{equation}The nonlinear time scale can be computed in our case as the spectral evolution time scale:
\begin{equation} \tau_{nl}=\frac{n_{k}}{\dfrac{\partial n_{k}}{\partial t}}, \end{equation}
with 
$\partial n_k/\partial t$ evaluated from our numerical data. Hereafter, we also use the terminology ‘spectral evolution time scale’ which is somewhat more illuminating than ‘nonlinear time scale’. In addition, we note that in many cases (e.g. Newell, Nazarenko & Biven Reference Newell, Nazarenko and Biven2006; Newell & Rumpf Reference Newell and Rumpf2011) the nonlinear time scale is taken as the kinetic time scale, which is not appropriate here since 
$\partial n_k/\partial t$ in (3.2) is evaluated by the dynamic equation instead of the WKE. The ratio of the spectral evolution and linear time scales is given by
\begin{equation} \rho=\frac{\tau_{nl}}{\tau_{lin}}. \end{equation}
The WKE is expected to be valid only when the time scales are well separated with 
$\rho \gg O(1)$.
To study the time-scale separation in our cases, we first mention that 
$\rho$ is, in general, a function of 
$k$, as shown in two typical cases in figure 4 for 
${\rm \Delta} \sigma =0.033$ and 
$1.461$. We see that 
$\rho$ generally increases with 
$k$ in a power-law form, which is related to our choice 
$\beta =0$ in (2.1) (which makes the nonlinear strength weaker for high 
$k$ than that from a positive 
$\beta$). Additionally, we see that as 
${\rm \Delta} \sigma$ increases from the former values to the latter, 
$\rho$ generally drops below 
$O(1)$, indicating the breakdown of the WKE in the latter case. We next seek to understand the relationship between the time-scale separation and the transition to the dynamic scaling of spectral evolution. To this end, we first note that the transition to dynamic scaling at all inertial-range wavenumbers (see figure 3) occurs at a similar forcing perturbation, indicating that the WKE breaks down for the overall spectral range instead of a particular wavenumber. Therefore, we evaluate 
$\rho$ at a representative wavenumber 
$k=100$, which leads to a relatively low value of 
$\rho$ in the inertial range as an assessment of the overall validity of wave turbulence theory. We plot in figure 5 
$\rho (k=100)$ as a function of 
$\varepsilon$, overlaid with the previous plot of 
$\partial n_{k}/\partial t$ for 
$k=500$ (as an example). We see that the transition to dynamic scaling occurs at 
$\rho \approx 4 \sim O(1)$, which indicates that the dynamic scaling range is consistent with the failure of the WKE due to the violation of time-scale separation. We also remark that if a different value of 
$k$ is chosen for the evaluation of 
$\rho$, we will end up with a slightly different transition value of 
$\rho$ which is at most approximately 10 as one can estimate from figure 4.
Figure 4. Ratio of time scales 
$\rho$ as a function of 
$k$ for cases with 
${\rm \Delta} \sigma =0.033$, 
$\varepsilon =0.022$ (green circles) and 
${\rm \Delta} \sigma =1.461$, 
$\varepsilon =0.124$ (red circles).
Figure 5. Ratio of the time scales 
$\partial n_k/\partial t(k=500)$ (left axis, blue triangles) and 
$\rho (k=100)$ (right axis, circles) as functions of 
$\varepsilon$. The region of transition from kinetic to dynamic scaling is circled, with the corresponding value of 
$\rho$ indicated (dash-dotted).
While the above analysis reveals the reduction of 
$\tau _{nl}$ as the underlying reason for the transition to dynamic scaling, the cause of this reduction remains unclear. In many previous studies, the reduction of 
$\tau _{nl}$ is attributed to the increase of nonlinearity level, which transits the dynamics into a strong turbulence regime (often associated with intermittency) where the WKE becomes irrelevant. However, we argue that this is not the major reason for the transition to dynamic scaling observed in our study, which is instead mainly triggered by the spectrum being too far from the stationary spectrum after a strong forcing perturbation is applied. To provide evidence for this argument, we repeat our analysis for three additional situations with higher values of 
$\varepsilon _0$. For each case, we evaluate a transition value 
$\varepsilon$ by the following procedure. We first define a dynamic scaling range as the longest range at a high nonlinearity level where a linear fit of 
$\varepsilon ^4$ has a correlation coefficient 
$R^2>0.998$. The transition value 
$\varepsilon$ is then defined as the leftmost point in the dynamic range. Figure 6 shows that the transition value of 
$\varepsilon$ increases with 
$\varepsilon _0$ following an approximately linear form of 
$\varepsilon \sim \varepsilon _0$, and that these transitions all correspond to 
$\rho (k=100) \sim O(1- 10)$, which is in agreement with the above analysis. This is consistent with the argument about the transition caused by a far-from-stationary spectrum induced by a strong forcing perturbation, since the deviation of the spectrum from the stationary state is approximately measured by the difference between 
$\varepsilon$ and 
$\varepsilon _0$, which does not change significantly in all cases. Moreover, figure 6 is in clear contradiction with a transition induced by a high nonlinearity level (or strong turbulence), in which case one would otherwise expect that the transition occurs at approximately the same 
$\varepsilon$ for different 
$\varepsilon _0$.
Figure 6. Values of 
$\varepsilon$ (left axis, squares) and 
$\rho (k=100)$ (right axis, circles) corresponding to transition to dynamic scaling for each base nonlinearity level 
$\varepsilon _{0}$, with a linear fitting (dotted line) and relation indicated.
The failure of the WKE due to the spectrum far from the stationary state can also be reasoned from a thought experiment. Given a spectrum, one can assume a priori that the WKE is valid, based on which 
$\partial n_{k}/\partial t$ can be computed. Such computed 
$\partial n_{k}/\partial t$ may a posteriori be found to violate the condition of 
$\rho \gg O(1)$ that invalidates the assumed WKE. A typical example of this violation has long been known in the Garrett–Munk spectrum of internal gravity waves, where the time-scale separation is not valid in the high-wavenumber/high-frequency portion of the spectrum (e.g. Holloway Reference Holloway1980; McComas & Müller Reference McComas and Müller1981; Lvov, Polzin & Yokoyama Reference Lvov, Polzin and Yokoyama2012; Eden, Pollmann & Olbers Reference Eden, Pollmann and Olbers2019; Wu & Pan Reference Wu and Pan2023). In this view, as the spectrum deviates from the stationary state, there exists a critical spectral form (even if the nonlinearity level is relatively low) that drives fast evolution beyond which the WKE fails.
3.2. Dominant interactions in modal growth
From § 3.1 we understand that the forcing perturbation (and the associated condensate) creates a deviation of the spectrum from the stationary state, which drives the subsequent spectral evolution. It is therefore reasonable to further argue that the inertial-range modal growth is dominated by direct interaction with the condensate peak, or at least non-local interactions with large-scale features. In this section, we verify this hypothesis using a new set-based interaction analysis. We remark that this is not a trivial task in the sense that these interactions cannot be analysed on the basis of the WKE, but rather must be analysed based on the dynamic equation (2.1) which is valid for all cases.
To start, we first define a set-based modal growth rate 
$\partial n_k/\partial t\rvert _{A}$ to be the evolution rate of 
$n_k$ due to the interaction of mode 
$k$ with 
$k_i\in A$ for 
$i=1,2,3$. By definition, if 
$A=\varLambda \equiv [-k_{max}, k_{max}]$ (i.e. the full spectral range with 
$k_{max}=1024$ after de-aliasing), then 
$\partial n(k)/\partial t\rvert _{A}=\partial n(k)/\partial t$. For 
$A\in \varLambda$, we can compute 
$\partial n(k)/\partial t\rvert _{A}$ by invoking (2.1) as
\begin{equation} \left.\frac{\partial n(k)}{\partial t}\right\rvert_{A} = \sum_{k_{1} \in A, k_{2}\in A, k_{3}\in A} 2\textrm{Im} \langle\hat{\psi}_{1}\hat{\psi}_{2}\hat{\psi}^\ast_{3}\hat{\psi}^\ast_{k}\rangle\delta^{12}_{3k}, \end{equation}
where 
$\delta ^{12}_{3k}\equiv \delta (k_1+k_2-k-k_3)$ and 
${\hat {\psi }}_k^*$ denotes the complex conjugate of 
$\psi _k$. The direct computation of (3.4) in spectral space is very expensive, with a computational complexity of 
$O(k_{max}^3)$ for each ensemble member of the group at each time instant (the ensemble and the time average are then used to compute the operator 
$\langle \cdots \rangle$). The computational cost can be significantly reduced if the evaluation of (3.4) can be performed in physical space, which turns out to be possible as we detail in Appendix A, with the result
\begin{equation} \frac{\partial n(k)}{\partial t}\biggr\rvert_{A}\equiv 2\textrm{Im} \langle \hat{\psi}^*_{k} \mathcal{F} [B_{A}(\psi(x))B_{A}(\psi(x))B_{A}(\psi^*(x))] \rangle, \end{equation}
where 
$\mathcal {F}[\cdots ]$ represents the Fourier transform and 
$B_A$ represents a Fourier-domain filter to select modes in 
$A$ for a quantity in the physical domain. The computation of (3.5) only requires a number of fast Fourier transforms which is much less expensive than that for (3.4).
In order to sort out the dominant interactions that lead to the growth of a mode 
$k_0$, we evaluate 
$\partial n(k_0)/\partial t\rvert _{A}$ by setting 
$A=[-k_{max}, -k_A]\cup [k_A, k_{max}]$ with 
$k_A$ varying in 
$[1,k_0]$. Therefore, for fixed 
$k_0$, 
$\partial n(k_0)/\partial t\rvert _{A}$ can be considered as a function of 
$k_A$. As 
$k_A$ decreases from 
$k_0$ to 1, 
$\partial n(k_0)/\partial t\rvert _{A}$ accounts for more non-local interactions with large scales (which is more important than interactions with small scales that are not particularly studied in the current setting). Figure 7 plots 
$\partial n(k_0)/\partial t\rvert _{A}$ (normalised by 
$\partial n(k_0)/\partial t$) as a function of 
$k_A$ for three selected values of 
$k_0$ in the inertial range, and for two cases of 
$\varepsilon =0.022$ and 
$0.124$ corresponding to regimes of kinetic and dynamic scaling, respectively. We first see that, for each case, the leftmost point (i.e. when 
$k_A=1$) recovers the total 
$\partial n(k_0)/\partial t$, which can be considered as a validation of our computation. For 
$k_A$ relatively close to 
$k_0$, local interactions are represented, with those for higher nonlinearity levels showing stronger fluctuations, i.e. stronger local interactions. Nevertheless, for all cases, there exists a wavenumber 
$k_c$ (marked in figure 7) at which the normalised 
$\partial n(k_0)/\partial t\rvert _{A}$ becomes close to zero or negative. This means that the overall local interactions in the range 
$[k_c, k_0]$ do not contribute much (or even contribute negatively) to the total 
$\partial n(k_0)/\partial t$. The major contribution to recover 
$\partial n(k_0)/\partial t$ occurs for a range immediately left of 
$k_c$, where the normalised 
$\partial n(k_0)/\partial t\rvert _{A}$ grows quickly from the minimum value to 1. This range features non-local interactions with scale separation of 
$k_c/k_0\sim O(0.1)$ (see table 1 for exact numbers) for all cases shown in the figure. The fact that 
$k_c$ increases with 
$k_0$ indicates that the large scales participating in the dominant interactions propagate to higher wavenumbers. For 
$k_0=150$, the non-local interactions mainly involve the forced condensate range up to about 
$k_c\sim O(10)$. For 
$k_0=500$ and 
$k_0=800$, the small-wavenumber range sufficiently filled with energy during the propagation of the perturbation becomes the dominant interacting modes.
Figure 7. Rate of spectral evolution 
$\partial n(k_0)/\partial t\rvert _{A}$, normalised by 
$\partial n(k_0)/\partial t$, as a function of 
$k_A$, for three wavenumbers (a,b) 
$k_0=150$, (c,d) 
$k_0=500$ and (e,f) 
$k_0=800$. The panels on the left (a,c,e) and on the right (b,d,f) are for cases with 
${\rm \Delta} \sigma =0.033$ and 
$1.461$ in the kinetic and dynamic scaling regimes, respectively. The positions of 
$k_c$ (where minimum value 
$\partial n(k_0)/\partial t\rvert _{A}$ is achieved) and 
$k_0$ are indicated in each panel.
Table 1. Values of 
$k_c$ and 
$k_c/k_0$ for all test cases.
In summary, the analysis here implies that, for all forcing perturbations considered, the inertial-range modal growth is mainly driven by the non-local interactions of the modes with large scales. The local interactions introduce some fluctuations that are stronger for higher-nonlinearity cases, but are sufficiently decayed when enlarging the spectral range around 
$k_0$ and before the dominant non-local interactions take place. We also remark that in many situations it is possible to derive a diffusion approximation (assuming the WKE is valid) for non-local interactions as done in McComas & Müller (Reference McComas and Müller1981) for three-wave systems and recently in Korotkevich et al. (Reference Korotkevich, Nazarenko, Pan and Shatah2023) for four-wave systems. However, for a (nonlinear) power-law dispersion relation, the diffusion approximation is only valid for dimension 
$d\ge 2$, since for 
$d=1$ there are no exact resonances when scale-separated interactions are considered, e.g. for four-wave systems we refer to two large and two small wavenumbers. Therefore, such diffusion approximation is not applicable to our system as far as exact resonances are considered.
4. Conclusion and discussions
In this paper, we numerically study the time scales of spectral evolution of nonlinear waves under perturbed forcing, in order to understand the previously identified ‘fast’ spectral evolution with the dynamic time scale contradicting the prediction from the WKE. Our study is conducted in the context of the 1-D MMT equation that mimics the wave turbulence of surface gravity waves free of other associated complexities. We show that both kinetic and dynamic time scales can become relevant for spectral evolution depending on the magnitude of the perturbed forcing applied to the stationary spectrum. The kinetic scaling occurs at weaker forcing perturbations, which transits to dynamic scaling as the forcing perturbation becomes strong enough, with the transition corresponding to the situation of the spectral evolution (or nonlinear) time scale dropping to the same level as the linear time scale. Such decrease in nonlinear time scale is mainly induced by the spectrum being far from the stationary state after the forcing perturbation is applied, instead of solely due to the usually considered cause of increasing nonlinearity level. Finally, through a new set-based interaction analysis, we find that the inertial-range spectral growth is dominated by its non-local interactions with the forced condensate or large scales sufficiently filled with energy throughout all nonlinearity levels.
Our study suggests that in assessing the validity of the WKE for transient spectra, the specific spectral forms need to be considered, in conjunction with the traditionally considered factors such as weak/strong turbulence and intermittency. For example, the mechanism understood in this paper may apply to the initial evolution of the wave spectrum, which is sometimes reported to be on the dynamic scale (Dyachenko et al. Reference Dyachenko, Trulsen, Krogstad and Socquet-Juglard2003) and sometimes precisely matches the prediction by the WKE (Zhu et al. Reference Zhu, Semisalov, Krstulovic and Nazarenko2022).
In addition, the current results in the MMT model have highlighted mechanisms that have physical relevance in surface gravity waves. We can now expect that in the event that wind forcing provides an adequately strong spectral perturbation and drives the spectrum far from equilibrium, the subsequent spectral evolution may follow the dynamic time scaling instead of prediction from WKE. The growth of spectral energy at small scales in this case should be expected to be governed by their non-local interactions with the forcing peak. Of course there exist many open questions regarding the case for gravity waves, e.g. How strong should the wind forcing and spectral perturbation be to realise dynamic scaling? Does wind direction relative to waves matter in this case? Does the three-wave process that generates bound modes play a role in the spectral evolution of gravity waves? We will consider a future study directly regarding surface gravity waves by applying the methodology developed in this paper, where some of the questions above will be answered.
Acknowledgements
The authors thank Professor M. Onorato for helpful discussions on this work during the Simons Collaboration on Wave Turbulence Annual Meeting in 2022.
Funding
The authors acknowledge the Simons Foundation for funding support for this work.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Fast computation for set-based interaction analysis
We start by reordering (3.4) (by interchanging allowed operations) as
\begin{equation} \frac{\partial n(k)}{\partial t}\biggr\rvert_{A} = 2\textrm{Im} \langle \hat{\psi}^*_{k} \sum_{k_{1} \in A, k_{2}\in A, k_{3}\in A} \hat{\psi}_{1}\hat{\psi}_{2}\hat{\psi}^*_{3} \delta^{12}_{3k} \rangle. \end{equation}The significant computational cost comes from the summation term on the right-hand side of (A1). To reduce the computational cost, we can consider the summation term as the Fourier transform of a physical-space quantity:
\begin{equation} \sum_{k_{1} \in A, k_{2}\in A, k_{3}\in A} \hat{\psi}_{1}\hat{\psi}_{2}\hat{\psi}^*_{3}\delta^{12}_{3k} = \mathcal{F} [B_{A}(\psi(x))B_{A}(\psi(x))B_{A}(\psi^*(x))], \end{equation}
where 
$\mathcal {F}[\cdots ]$ represents the Fourier transform and 
$B_A$ represents a Fourier-domain filter to select modes in 
$A$ for a physical-domain function 
$\psi (x)$. Equation (A2) can be understood from the convolution theorem for three functions, i.e. spectral-domain convolution is equal to the physical-domain multiplication. One can also directly confirm that this is true by starting from the right-hand side and expressing each 
$B_{A}(\psi (x))$ as the summation of Fourier modes:
\begin{equation} \frac{1}{2{\rm \pi}}\int \exp({-{\rm i}kx}) \sum_{k_1\in A} \hat{\psi}(k_1) \exp({{\rm i}k_1x}) \sum_{k_2\in A} \hat{\psi}(k_2) \exp({{\rm i}k_2x}) \sum_{k_3\in A} \hat{\psi}^*(k_3) \exp({-{\rm i}k_3x})\,{{\rm d}\kern0.06em x}, \end{equation}which can be reorganised as
\begin{equation} \frac{1}{2{\rm \pi}}\int \sum_{k_{1} \in A, k_{2}\in A, k_{3}\in A} \hat{\psi}(k_1) \hat{\psi}(k_2) \hat{\psi}^*(k_3) \exp({{\rm i}(k_1+k_2-k-k_3)x})\, {{\rm d}\kern0.06em x}. \end{equation}
We see that (A4) is the same as the left-hand side of (A1) after considering 
$({1}/{2{\rm \pi} })\int \exp ({\textrm {i}(k_1+k_2-k-k_3)x}) = \delta (k_1+k_2-k-k_3)$.
Appendix B. Results from cases with focusing nonlinearity
In this appendix, we summarise results from cases with focusing nonlinearity, i.e. 
$\lambda =-1$ in (2.1). All other parameters in the study are kept consistent with those in the defocusing cases reported in the main paper. Figures 8 and 9 show respectively the scaling of 
$\partial n_k/\partial t$ with 
$\varepsilon$ and the study regarding 
$\rho$, as counterparts of figures 3 and 4 in the main paper. It is clear from these figures that the main conclusion made for the defocusing case also applies to the focusing case.
Figure 8. See caption for figure 3, but for the focusing case.
Figure 9. See caption for figure 4, but for the focusing case.
Open access
