In many-body physics, a central topic is the thermalization of systems far from equilibrium under Hamiltonian dynamics, especially when the Hamiltonian is of the form H = H₀ + gV. Typically, H₀ is integrable, while gV introduces non-integrability. While there is a consensus on the occurrence of thermalization [1–3], the route to thermalization is layered with complexity. The foundational work of Fermi, Pasta, Ulam, and Tsingou (FPUT) in the 1950s laid the groundwork for this topic [4]. Their insights paved the way for a cascade of studies that have progressively deepened our grasp [5–12]. Recently, multi-wave resonances have been spotlighted as the driver for energy redistribution, culminating in thermalization [13–16]. This process, encapsulated by the Zakharov equation [17, 18], suggests that systems reach thermalization on a timescale proportional to 1/g² [19–22]. In the quantum realm, a parallel story unfolds, with the quantum Boltzmann equation superseding its classical counterpart [23–25]. Theoretical success has been achieved when the perturbation is weak, i.e., g ≪ 1. Such a perturbation ensures that the system quickly settles to a quasi-stationary state encoded via a generalized Gibbs ensemble (GGE) [26] and relaxes to the final thermal state on a much larger timescale ∝1/g² (See figure 1). This separation of timescale has been termed prethermalization, a subject that has garnered significant attention [27–32].

However, while weak interactions have been extensively studied, the domain of strong interactions remains a frontier, teeming with unanswered questions. Preliminary research, such as the work by Lepri [6], has hinted at the profound implications of strong interactions, but the underlying mechanisms remain elusive. What role do multi-wave resonances assume in this regime? Does prethermalization, a hallmark of weakly interacting systems, still manifest? Exploring this realm is crucial for understanding many-body systems and has implications for its applications, including light propagation in multi-mode fibers [33, 34] and energy cascades in wave-turbulence systems [35].

In this paper, we delve deep into the realm of strong interactions, a territory that holds the promise of reshaping our understanding of thermalization. We introduce a new framework for studying thermalization in the system of interacting waves. Distinguishing between trivial (no momentum exchange) and nontrivial interactions (significant energy redistribution), we transition to a model of interacting dressed waves with weakened interactions. Exemplified by FPUT-like models, we show that prethermalization survives the strong interaction and that thermalization is still dominated by resonances among waves, which are now effectively dressed. This dominance, coded in the Zakharov equation, leads to a double scaling in the thermalization time, further supported by our numerical results.

Our framework addresses systems of interacting waves or quasi-particles. Utilizing the vector a to denote the collection of variables {a₁, a₂, …} and their complex conjugates $\{{a}_{1}^{* },{a}_{2}^{* },\ldots \}$, the governing Hamiltonian is given bywhere ω_k represents the frequency spectrum and k signifies momentum modes. The unperturbed Hamiltonian, typically integrable, is denoted by H₀. The term responsible for breaking integrability arises from p-wave interactions (p > 2):with its magnitude determined by g. For simplicity, we employ shorthand notations: j = (j₁, ⋯ ,j_q), and l = (l₁, ⋯ ,l_r), leading to ${V}_{{\boldsymbol{j}}}^{{\boldsymbol{l}}}={V}_{{j}_{1}\ldots {j}_{q}}^{{l}_{1}\ldots {l}_{r}}$. The interaction coefficient, owing to the Hamiltonian's nature, exhibits the symmetry ${V}_{{\boldsymbol{j}}}^{{\boldsymbol{l}}}={{V}_{{\boldsymbol{l}}}^{{\boldsymbol{j}}}}^{* }$. For translation-invariant interactions, ${V}_{{\boldsymbol{j}}}^{{\boldsymbol{l}}}$ is nonzero only when the momentum conservation condition ${k}_{{j}_{1}}+\cdots +{k}_{{j}_{q}}={k}_{{l}_{1}}+\cdots +{k}_{{l}_{r}}$ is satisfied. Drawing inspiration from quantum mechanics, we interpret a_k as incoming waves and ${a}_{k}^{* }$ as outgoing waves. The term ${a}_{{j}_{1}}\cdots {a}_{{j}_{q}}{a}_{{l}_{1}}^{* }\cdots {a}_{{l}_{r}}^{* }$ can be visualized as a process where r outgoing waves arise from the interaction among q incoming waves.

Given an initial state f_I, the isolated system evolves for t > 0 as: f(t) = e^−Ltf_I. Here, L is the Liouville operator, defined by the Poisson bracket as L = {H,} . Over time, the system is expected to achieve equilibrium, represented by a Gibbs ensemble (GE) at temperature θ:where $\mathrm{Tr}[f]$ is the partition function.

To describe the post-equilibration expectation values of observables when H is nonintegrable, one typically substitutes the GE (3) with a Gaussian distribution [36], with mean vector α = ⟨a⟩_GE and covariance matrix Θ = ⟨(a^† − α^†)(a − α)⟩_GE. Here, the superscript ‘†' represents the conjugate transpose, and the ‘GE' subscript indicates an average taken under GE (3). The formulation (4) is rooted in mean-field theory [36]. Without loss of generality, we will focus on cases where α = 0.

The Hermiticity of the matrix Θ implies the existence of a unitary matrix U such that ${[{U}^{\dagger }{\rm{\Theta }}U]}_{{kl}}=\theta {\delta }_{{kl}}/{\tilde{\omega }}_{k}$. By transitioning to a new set of coordinates, $\tilde{{\boldsymbol{a}}}=U{\boldsymbol{a}}$, the distribution (4) can be expressed as:This distribution represents the so-called GGE of an integrable Hamiltonian system:which characterizes a system of noninteracting dressed waves.

The general equipartition theorem says that ⟨a_k∂H/∂a_l⟩_GE = θδ_kl and $\langle {a}_{k}\partial H/\partial {a}_{l}^{* }{\rangle }_{\mathrm{GE}}=0$ [37]. From this, we deduce For a system of noninteracting waves (i.e., g = 0), the waves remain free, and Θ is diagonal. However, the introduction of nonlinear interaction causes waves to become correlated, making Θ generally non-diagonal. If one expands GE as $\exp (-H/\theta )=\exp (-{H}_{0}/\theta )[1-g/\theta {V}^{(p)}({\boldsymbol{a}})+\cdots ]$, it is ready to conclude that not all interactions contribute to Θ. This insight enables us to categorize the interactions V^(p)(a) into two groups: (i) ${V}_{0}^{(p)}$: these interactions involve pairwise equal momenta for incoming and outgoing waves, e.g., ${a}_{k}\cdots {a}_{l}{a}_{k}^{* }\cdots {a}_{l}^{* }$. In other words, trivial interactions occur without any momentum exchange. (ii) ${V}_{\mathrm{non}}^{(p)}$: these interactions involve momentum exchange. In summary, ${V}^{(p)}\,={V}_{0}^{(p)}+{V}_{\mathrm{non}}^{(p)}$. The contribution of ${V}_{0}^{(p)}$ to the covariance Θ is clearly evident, while ${V}_{\mathrm{non}}^{(p)}$ does not exhibit such an influence. As discussed in more detail in [38], interactions that facilitate energy redistribution are consistently classified as nontrivial.

The explicit expression for ${H}_{0}(\tilde{{\boldsymbol{a}}})+{V}_{0}^{(p)}(\tilde{{\boldsymbol{a}}})$ depends on the unitary transformation ${\boldsymbol{a}}={U}^{\dagger }\tilde{{\boldsymbol{a}}}$. While ${\tilde{H}}_{0}(\tilde{{\boldsymbol{a}}})\,={\sum }_{k}{\tilde{\omega }}_{k}{\tilde{a}}_{k}{\tilde{a}}_{k}^{* }$ and ${H}_{0}(\tilde{{\boldsymbol{a}}})+{V}_{0}^{(p)}(\tilde{{\boldsymbol{a}}})$ are formulated differently, our analysis reveals that ${\tilde{H}}_{0}$ contains the statistical information embedded in H₀ and ${{gV}}_{0}^{(p)}$. By reorganizing the original Hamiltonian (1), we derive a Hamiltonian that is statistically equivalent to it:This effective Hamiltonian shares the same form as the original Hamiltonian H. The primary distinction between them lies in that $\tilde{H}$ excludes the trivial interactions. Statistically, $g{\tilde{V}}_{\mathrm{non}}^{(p)}$ does not significantly influence the energy landscape defined by ${\tilde{H}}_{0}$. The exclusion of trivial interactions results in $g\langle {V}_{\mathrm{non}}^{(p)}\rangle \approx 0$. Consequently, ${{gV}}_{\mathrm{non}}^{(p)}$ primarily serves as a medium for energy transfer and, on average, does not contribute significantly to the total energy. In this sense, we map the original system (1) into a more weakly interacting one. For such a system (8), at the large timescales ($t\,\gg {\tau }_{0}\sim 2\pi /{\tilde{\omega }}_{k}$), the phase-space distribution $\tilde{f}(t)$ can be approximated by a GGE:Due to statistical equivalence, we expect that the GGE (9) can further replace the phase-space distribution f(t) of the original system. If the system thermalizes, one expects that θ_k(t) = θ and the distribution (5) results. This phenomenon is often referred to as prethermalization. We outline our understanding of the prethermalization framework in figure 1. In systems exposed to abrupt perturbations, statistically uncorrelated quasi-particles swiftly emerge from the system's complex dynamics. Initially, these particles exhibit a spectrum of temperatures, denoted as θ_k(t). However, through subsequent collisions, they converge to a consistent temperature over large timescales. This behavior suggests that a Boltzmann-like equation might aptly describe this evolution, a proposition we will explore in depth in the subsequent sections.

The transformation $H\to \tilde{H}$ offers an alternative approach to discern the dynamics of observables O. Next, we provide the second-order approximation to ∂_t⟨O⟩. To achieve this, we decompose the Liouville operator of the equivalent system (8):Making necessary approximations [39], we obtainwhere ${\tilde{L}}_{V}(t)={{\rm{e}}}^{-{\tilde{L}}_{0}t}{\tilde{L}}_{V}{{\rm{e}}}^{{\tilde{L}}_{0}t}$. We focus on situations where the initial distribution is ${\tilde{f}}_{{\rm{I}}}$ such that ${\tilde{L}}_{0}{\tilde{f}}_{{\rm{I}}}=0$. Even if the initial state is not set up this way, prethermalization will rapidly steer it to a state that fulfills this condition [1]. A direct calculation yields the time-dependent average of observable O:where $\tilde{f}(t)$ has been replaced with the GGE (9) when averaging.

To elucidate our findings and align with prevalent numerical experiments, we employ the FPUT-like model. The Hamiltonian is expressed as:where ${\dot{u}}_{j}$ and u_j represent the velocity and displacement of the jth particle, respectively. In this equation, the quadratic terms form the unperturbed Hamiltonian H₀, while the nonlinear terms break the integrability. The order of nonlinearity p must be even to ensure the stability in strongly interacting regimes.

Transformations (see appendix A) recast the Hamiltonian into the form (1). The interaction strength correlates with energy density ϵ = E/N and nonlinearity strength λ, specifically as g = λϵ^(p−2)/2. For the FPUT-like model, nonlinearity induces correlations between waves a_k and a_N−k [38, 40]. The transformation that diagonalizes matrix Θ is (see appendix B.1):with η termed the dressing/renormalization factor, whose asymptotic behavior outlined in the appendix B.2. The statistically effective Hamiltonian iswhere ${\tilde{\omega }}_{k}=\eta {\omega }_{k}$. The interaction ${V}_{\mathrm{non}}^{(p)}$ excludes trivial interactions, effectively reducing nonlinear interactions (η ≥ 1) [38, 40, 41].

In the context of thermalization, our primary focus is on the dynamics of the wave action spectral density, denoted as ${\tilde{n}}_{k}=\langle {\tilde{a}}_{k}^{* }{\tilde{a}}_{k}\rangle $. By substituting $O={\tilde{a}}_{k}^{* }{\tilde{a}}_{k}$ into equation (12), we derive the desired Boltzmann-like equation, which is frequently termed the Zakharov equation:The explicit expressions of coefficients ζ_k and γ_k are detailed in the appendix C. The p-wave resonance conditions are contained in both ζ_k and γ_k, i.e., k₁ ± ⋯ ± k_p = 0 and ${\tilde{\omega }}_{{k}_{1}}\pm \cdots \pm {\tilde{\omega }}_{{k}_{p}}=0$. Only when both ζ_k and γ_k are nonvanishing, does the kinetic equation (16) suffice. Otherwise, canonical transformations should be introduced to explore higher-order multi-wave resonances [17, 18]. The resonance conditions imply that energy is transferred among resonant waves [19–22]. When a wave engages in multiple resonances, it links these resonances, allowing energy to cascade through a larger network of waves. This phenomenon, where waves/resonances are ‘connected' through wave sharing, facilitates effective energy redistribution across the entire system. Chirikov's resonance-overlap criterion offers an alternative theoretical framework for understanding the connectivity of resonances [42]. It posits that chaos and extensive energy redistribution can occur when two resonances significantly overlap. While our model primarily relies on the simpler mechanism of direct wave sharing to connect resonances, Chirikov's criterion provides a potential explanation for more complex interactions where wave sharing alone might be insufficient.

The kinetic equation (16) obtained after the averaging procedure, however, describes an irreversible evolution towards thermodynamic equilibrium. The mathematical statement of irreversibility is the theorem of entropy growth which is similar to Boltzmann's H-theorem for gas kinetics. Indeed, let us consider the time dependence of the entropywhere we have used the relation ${\tilde{n}}_{k}=\langle {\theta }_{k}(t)\rangle /{\tilde{\omega }}_{k}$. After some straightforward calculations, they readily obtain dS/dt ≥ 0. Thus, the entropy of a closed wave system can only increase. The thermodynamic equilibrium corresponds to the maximum of entropy given some constraints. Indeed, it is easy to see that the energy ${ \mathcal E }={\sum }_{k}{\tilde{\omega }}_{k}{\tilde{a}}_{k}{\tilde{a}}_{k}^{\ast }$ is an integral of motion in equation (16). Making use of the method of Lagrange multipliers, we obtainThe parameter β is given by the temperature of the system, i.e., β = 1/θ, hence we have the Rayleigh–Jeans distribution ${\tilde{n}}_{k}=\displaystyle \frac{\theta }{{\tilde{\omega }}_{k}}$, which corresponds to the energy equipartition, i.e., ${\tilde{\omega }}_{k}{\tilde{n}}_{k}=\theta $.

Our findings corroborate that equation (16) is consistent with thermalization. Multi-wave resonances are instrumental in energy redistribution and eventual thermalization. From equation (16), the thermalization time is inferred to scale as [21, 22]:The validity of this result hinges on the divergence of linear and nonlinear timescales, specifically when $2\pi /{\tilde{\omega }}_{k}\ll 1/{\gamma }_{k}$ [17, 18]. Our result improves the narrative of thermalization. Specifically, it unveils that it is the resonances among dressed waves ${\tilde{a}}_{k}$ that govern the thermalization. In the weakly interacting regime (WIR) where η ∼ 1, τ_th ∝ g⁻², a conclusion supported by most recent works [13–15, 19–22]. Conversely, in the strongly interacting regime (SIR) with η ∼ g^1/p, the thermalization timescales as τ_th ∝ g^−1/p. Summarily,Concerning the SIR, the asymptotic behavior of the thermalization time can also be derived using a simple scaling argument in the high-energy limit. Specifically, in this limit, one can normalize the total energy to unity through the following rescalings: u_j → (Nϵ)^1/pλ^−1/pu_j, t → λ^−1/p(Nϵ)^1/p−1/2t and${\dot{u}}_{j}\to {(N\varepsilon )}^{1/2}{\dot{u}}_{j}$.

In our subsequent analysis, we numerically validate our theoretical insights. For the FPUT-like system, the thermalized state is described by equation (5), leading to energy equipartition, i.e., ${\tilde{\omega }}_{k}\langle {\tilde{a}}_{k}^{* }{\tilde{a}}_{l}\rangle =\theta {\delta }_{{kl}}$. To gauge thermalization, we employ the effective degrees of freedom introduced in [11], defined as d_eff(t) = 2ξ(t)e^s(t)/N, where $s(t)=-{\sum }_{k=N/2}^{N}{w}_{k}(t)\mathrm{ln}{w}_{k}(t)$ is the spectral entropy, with ${w}_{k}={\bar{e}}_{k}(t)/{\sum }_{l=N/2}^{N}{\bar{e}}_{l}(t)$, $\xi (t)\,=2{\sum }_{k=N/2}^{N}{\bar{e}}_{k}(t)/{\sum }_{k=1}^{N}{\bar{e}}_{k}(t)$ and ${\bar{e}}_{k}(t)\,=1/(T-\mu T){\int }_{\mu T}^{T}{e}_{k}(t){\rm{d}}t$. Here, e_k(t) could be the energy of either a_k or ${\tilde{a}}_{k}$, and parameter μ controls the size of the time window for averaging and is fixed at μ = 2/3 in computation. In the upcoming analysis, we contrast the effective degrees of freedom derived from a_k and ${\tilde{a}}_{k}$, denoted as d_eff(a, t) and ${d}_{\mathrm{eff}}(\tilde{{\boldsymbol{a}}},t)$ respectively. We anticipate ${d}_{\mathrm{eff}}(\tilde{{\boldsymbol{a}}},t)=1$ when the system thermalizes.

For our simulations, we utilize the eighth-order Yoshida algorithm [43], adjusting the time step based on energy density to maintain an absolute error below 10⁻⁴. The initial energy is uniformly distributed among 10% waves of the lowest-frequency. Observables are averaged over 128 random initial states to mitigate fluctuations. We set λ = 1 and explore the dependence of τ_th on ϵ [See equations (20)].

Before moving on, we first clarify a question: can waves a_k be used to compute the levels of thermalization? Our findings affirmatively answer this. While equipartition within dressed waves ${\tilde{a}}_{k}$ consistently indicates thermalization, the relation between energies of a_k and ${\tilde{a}}_{k}$, namely, ${\omega }_{k}\langle {a}_{k}^{* }{a}_{k}{\rangle }_{\mathrm{GE}}/{\tilde{\omega }}_{k}\langle {\tilde{a}}_{k}^{* }{\tilde{a}}_{k}{\rangle }_{\mathrm{GE}}=({\eta }^{2}+1)/(2{\eta }^{2})$ under GGE ansatz (9), suggests that their thermalization levels are congruent. We validate this with numerical experiments on system ${{ \mathcal H }}_{4}$ with N = 1024. Figure 2 displays the effective degrees of freedom d_eff(t) for two energy densities, representing WIR and SIR. The figure reveals that thermalization levels among a (red circle) or $\tilde{{\boldsymbol{a}}}$ (blue cross) are nearly identical, corroborating our hypothesis (9).

We now delve into the characteristics of d_eff(t). Numerical experiments are conducted on the system ${{ \mathcal H }}_{4}$ with N = 1024. In figure 3, we depict d_eff(t) against time scaled by energy density (whose value is already listed in the figure) ϵ^αt, where α = 2 in the WIR (left) and α = 1/4 in the SIR (right). Notably, all d_eff(t) increase monotonically from 0 to 1, and data collapse is evident across different energy densities in both regimes, suggesting distinct thermalization timescales for WIR and SIR. We define the thermalization time when d_eff(t) reaches a threshold d_c. While the specific value of d_c does not influence the scales, it should be sufficiently large to capture significant thermalization in the SIR. In our computation, we set d_c = 0.9.

Figure 4 presents the thermalization time τ_th against energy density ϵ for systems ${{ \mathcal H }}_{4}$ (left) and ${{ \mathcal H }}_{6}$ (right). Theoretical predictions for WIR and SIR are represented by solid and dashed lines, respectively. The predictions align well with the data, but some deviations in the WIR are attributed to the finite-size effect, as reported in previous studies [19, 21]. Larger system sizes yield better agreement with the prediction τ_th ∝ ϵ⁻². To illustrate this, the computation with a larger size, i.e., N = 8192, is performed. Notably, strong interactions negate the finite-size effect due to significant frequency broadening. The results for system ${{ \mathcal H }}_{6}$ also align with predictions, i.e., τ_th ∝ ϵ⁻⁴ in the WIR and τ_th ∝ ϵ^−1/3 in the SIR. Previous observations in smaller ${{ \mathcal H }}_{4}$ [15] noted the double scaling of thermalization time. The study proposed a crossover in timescale from τ_th ∝ ϵ⁻⁴ to τ_th ∝ ϵ⁻¹ with increasing interaction. In smaller systems, six-wave resonances primarily drive thermalization, resulting in τ_th ∝ ϵ⁻⁴. However, the energy density at which τ_th ∝ ϵ⁻¹ in [15] is around 1, distant from the SIR. As interactions intensify, it is anticipated that the thermalization time will scale as τ_th ∝ ϵ^−1/4. By measuring the broadening of the frequencies, [15] argued that Chirikov's criterion could be used to identify the transition region. It is noteworthy that this work applied Chirikov's criterion to two different waves rather than to two nonlinear resonances, as would typically be expected in such analyses. Despite this methodological deviation, their findings are consistent with ours. Before the transition region, overlaps are rare, and waves a_k remain nearly free, i.e., η ∼ 1. However, as the energy density increases and crosses this regime, overlaps become frequent and significantly enhance renormalization. Consequently, a different timescale emerges. To verify this, we have plotted the numerical and analytical results for the normalization factor in appendix B.2.

In this work, we introduced a novel approach to address non-equilibrium dynamics in systems of interacting waves. Central to our method is the transformation of an interacting system into a statistically equivalent counterpart, albeit with reduced interactions. This transformation offers a fresh perspective in discerning the dynamics of observables inherent to the original system. Through this lens, we discovered the occurrence of prethermalization in both the weakly interacting regime and the strongly interacting regime. Furthermore, our analysis underscores the pivotal role of resonances among dressed waves in steering the system towards thermalization. By leveraging second-order perturbation theory on the statistically equivalent system, we derived the kinetic equation, which elucidates the irreversible route to equilibrium. Notably, this equation revealed a distinct scaling behavior for the thermalization time in the SIR.

However, many questions still need to be addressed. In our study, we consider only the typical initial states, which rapidly settle into a quasi-stationary state encoded by GGE. What would have happened if this was not the case [44]? Additionally, the moderately nonlinear regime, where g⁻² ∼ τ₀, remains largely unexplored; in this regime, hydrodynamic effects may become significant.

Our methodology holds promise as a tool to investigate many-body localization in disordered systems. The dual effects of trivial interactions—amplifying disorder by dressing the Anderson localized modes and concurrently attenuating the interactions statistically—might pave the way for stable localized quasi-integrals of motion. As a concluding remark, it is worth emphasizing that our approach can seamlessly transition to a quantum framework, opening avenues for further exploration in quantum systems.

We begin with the widely recognized Fourier transformation:It is important to note that ${Q}_{k}^{* }={Q}_{N-k}$ and ${P}_{k}^{* }={P}_{N-k}$. After this transformation, the Hamiltonian (13) becomeswithTo investigate wave interactions in such systems, canonical complex normal variables a_k are usually introduced aswhere ${\omega }_{k}=2\sin (k\pi /N)$ is the dispersion relation. Using a_k, the FPUT-like model's Hamiltonian is represented aswhere symbol j is shorthand for a set of indexes j₁, ⋯ ,j_q, and similar definitions for l. The term ${a}_{{j}_{1}}\cdots {a}_{{j}_{q}}{a}_{{l}_{1}}^{* }\cdots {a}_{{l}_{r}}^{* }$, weighted by ${V}_{{\boldsymbol{j}}}^{{\boldsymbol{l}}}$, can be interpreted as a process, denoted by (q → r), in which r outgoing waves are ‘created' as a result of interacting among q incoming waves. Here, ${V}_{{\boldsymbol{j}}}^{{\boldsymbol{l}}}$ reads explicitlywith ${{\rm{\Delta }}}_{{\boldsymbol{j}}}^{{\boldsymbol{l}}}={{\rm{\Delta }}}_{{j}_{1}+\cdots +{j}_{q}-{l}_{1}-\cdots -{l}_{r}}^{(N)}$. For convenience, we rescale Hamiltonian (A5) by energy density ϵ = H/N: ${a}_{k}\to \sqrt{N\varepsilon }{a}_{k}$ and H → NϵH. The result readswhich is exactly the Hamiltonian discussed in our main text [See equations (1) and (2)]. Here, $g=\lambda {\varepsilon }^{\tfrac{p-2}{2}}$ represents the interaction strength.

For fixed boundary condition where u_N+1 = u₀ = 0, the transformation from $({\boldsymbol{u}},\dot{{\boldsymbol{u}}})$ to (Q, P) readsLet ${ \mathcal N }=2N+2$ and use the identity $\sin \theta ={\rm{i}}({{\rm{e}}}^{{\rm{i}}\theta }-{{\rm{e}}}^{{\rm{i}}\theta })/2$ to rewrite the above transformation aswhere ${{ \mathcal Q }}_{k}=-{{ \mathcal Q }}_{{ \mathcal N }-k}=-{{\rm{i}}{Q}}_{k}$ and ${{ \mathcal P }}_{k}=-{{ \mathcal P }}_{{ \mathcal N }-k}=-{{\rm{i}}{P}}_{k}$, therefore ${{ \mathcal Q }}_{k}={{ \mathcal Q }}_{{ \mathcal N }-k}^{* }$ and ${{ \mathcal P }}_{k}={{ \mathcal P }}_{{ \mathcal N }-k}^{* }$. Again, we introduce the canonical complex normal variableswhere ${\omega }_{k}=\sin (k\pi /{ \mathcal N })$. Transforming with equation (A8) and then rescaling also lead to the equation (A6), merely by replacing N with ${ \mathcal N }$. This result implies that a (2N + 2)-particle FPUT-like model with fixed boundary condition is equivalent to an N-particle one with periodic boundary condition. For this reason and in order to suppress the finite-size effect, we adopt fixed boundary condition in numerical experiments.

Substituting equation (14) into equation (A6), or substituting equation (B3) into equation (A2) effectively gives the Hamiltonian that we expected, i.e.,where the prime denotes the summation that neglects the trivial interactions, which are collected into ${V}_{\mathrm{non}}^{(p)}$ in our main text [See equation (15)]. For clarity, next we are going to focus on the very effect of reciprocal interactions of the types (q ⇌ r) in equation (A6), because, as we will immediately see, interactions of the type (q ⇌ r) and $(q^{\prime} \rightleftharpoons r^{\prime} )$ are statistically independent under the GGE ansatz (5) if $q\ne q^{\prime} $. Specifically, the Hamiltonian of interest readswhich differs from equation (C2) in that the sum over q (or r) is omitted.

For now, we could apply equation (12) into the observable $O={\tilde{a}}_{k}{\tilde{a}}_{k}$ which commutes with ${\tilde{H}}_{0}(\tilde{{\boldsymbol{a}}})$. Obviously, $\langle {\tilde{L}}_{0}{\tilde{I}}_{k}\rangle =0$. The first-order approximation readswhere ${J}_{{j}_{1},{j}_{2},\ldots ,{j}_{q}}^{{l}_{1},{l}_{2},\ldots ,{l}_{r}}\equiv \langle {\tilde{a}}_{{j}_{1}}{\tilde{a}}_{{j}_{2}}\cdots {\tilde{a}}_{{j}_{q}}{\tilde{a}}_{{l}_{1}}^{* }{\tilde{a}}_{{l}_{2}}^{* }\cdots {\tilde{a}}_{{l}_{r}}^{* }{\rangle }_{{\tilde{f}}_{0}}$ is the p-th order correlation function. We restrict our attention to only those situations where the initial phases are random. Therefore, ${J}_{{j}_{1},{j}_{2},\ldots ,{j}_{q}}^{{l}_{1},{l}_{2},\ldots ,{l}_{r}}=0$ if q or p − q is odd, and ${J}_{{j}_{1},{j}_{2},\ldots ,{j}_{q}}^{{l}_{1},{l}_{2},\ldots ,{l}_{r}}={{J}_{{l}_{1},{l}_{2},\ldots ,{l}_{r}}^{{j}_{1},{j}_{2},\ldots ,{j}_{q}}}^{* }$ if q and p are both even. Keeping in mind that ${W}_{{j}_{1},{j}_{2},\ldots ,{j}_{q}}^{k,{l}_{2},\ldots ,{l}_{r}}={{W}_{k,{l}_{2},\ldots ,{l}_{r}}^{{j}_{1},{j}_{2},\ldots ,{j}_{q}}}^{* }$, these properties indicate that up to the first-order perturbation we have ${\tilde{n}}_{k}=0$. We now turn to the second-order perturbation, i.e., ${\int }_{0}^{t}\langle {\tilde{L}}_{V}{\tilde{L}}_{V}(t-s)O{\rangle }_{\mathrm{GGE}}{\rm{d}}s$. After some straightforward and tedious calculations, we havewhereIn the calculation, we have used Wick's rule, which works for the GGE ansatz (5), and identitiesThe equation (C4) can be divided into two parts according to whether ${\tilde{n}}_{k}$ is included, leading the result (16). The result (16) or (C4) can also be obtained in the framework of wave turbulence approach [17, 18].