2. Review of normal form for $\beta$-FPUT model

We first briefly review the formulation of the $\beta$-FPUT model and its associated normal form transformation following [Reference Bustamante, Hutchinson, Lvov and Onorato3]. The $\beta$-FPUT model consists of $N$ identical masses $m$ connected by non-linear springs with the elastic force $F=-\kappa\Delta q+\beta\Delta q^3$, where $\kappa,\beta\neq 0$ are spring elastic constants. The equation of motion for a $\beta$-FPUT chain can be expressed as (e.g., [Reference Bustamante, Hutchinson, Lvov and Onorato3, Reference Fermi, Pasta, Ulam and Tsingou13]):

(2.1)\begin{equation} m\ddot{q}_j=\kappa\left[(q_{j+1}-q_j)-(q_j-q_{j-1})\right]+\beta\left[(q_{j+1}-q_j)^3-(q_j-q_{j-1})^3\right], \end{equation}

where $j=0,1,...,N-1$ and $q_j(t)$ represents the displacement from the equilibrium position of the mass $j$ with periodic boundary conditions $q_0=q_N$.

We consider a discrete Fourier transform

(2.2)\begin{equation} Q_k=\frac{1}{N}\sum_{j=0}^{N-1}q_je^{-i2\pi k_j/N},\mbox{ }P_k=\frac{1}{N}\sum_{j=0}^{N-1}p_je^{-i2\pi k_j/N}, \end{equation}

with $p_j=m\dot{q}_j$ and introduce the normal modes $a_k(t)$ as

(2.3)\begin{equation} a_k=\frac{\sqrt{2}}{2}\left[(m\omega_k)^{\frac{1}{2}}Q_k+i(m\omega_k)^{-\frac{1}{2}}P_k\right], \mbox{ for }k=1,...,N-1. \end{equation}

The evolution equation (2.1) then becomes

(2.4)\begin{align} i\dot{a}_{k_1}=&\omega_{k_1}{a}_{k_1}+\frac{\beta}{3}\sum_{k_2,k_3,k_4}\bigg[T_{1,2,3,4}a_{k_2}a_{k_3}a_{k_4}\delta(k_1-k_2-k_3-k_4)\\ notag &+3T_{1,-2,3,4}a_{k_2}^*a_{k_3}a_{k_4}\delta(k_1+k_2-k_3-k_4)+3T_{1,-2,3,4}a_{k_2}a_{k_3}^*a_{k_4}^*\delta(k_1-k_2+k_3+k_4)\\ notag &+T_{1,2,3,4}a_{k_2}^*a_{k_3}^*a_{k_4}^*\delta(k_1+k_2+k_3+k_4)\bigg], \end{align}

where the frequency and interaction coefficients read

(2.5)\begin{equation} \omega_{k}=\omega(k)=2\sqrt{\frac{\kappa}{m}}\left|\sin\left(\frac{\pi k}{N}\right)\right|, \end{equation}

(2.6)\begin{equation} T_{1,2,3,4}=- \frac{3}{4\kappa^2}\iota( k_2+k_3+k_4)\iota(k_2)\iota(k_3)\iota(k_4) \prod_{i=1}^4 \sqrt{\omega_{k_i}}, \end{equation}

with $\iota(x):=\text{sgn} \sin \left(\frac{\pi x}{N}\right)$ a $2N$-periodic function. $\delta(\cdot)$ denotes the Kronecker delta with mod $N$ applied to its argument. Hereafter the mod $N$ condition is applied to all wavenumber conditions but sometimes not written explicitly for simplicity.

Among the four terms in (2.4), only the second one with $\delta(k_1+k_2-k_3-k_4)$ leads to resonant interactions. For the other three non-resonant terms, the corresponding frequency condition cannot be satisfied due to the positivity and strict subadditivity of $\omega(k)$, as well explained in [Reference Bustamante, Hutchinson, Lvov and Onorato3]. Therefore, we can remove these non-resonant terms at the current order via a normal form transformation written as

(2.7)\begin{align} a_{k_1} =& b_{k_1}+\frac{\beta}{3}\sum_{k_2,k_3,k_4}\bigg[A^{(1)}_{1,2,3,4}b_{k_2}b_{k_3}b_{k_4}\delta(k_1-k_2-k_3-k_4)\\ notag &+3A^{(2)}_{1,2,3,4}b_{k_2}b_{k_3}^*b_{k_4}^*\delta(k_1-k_2+k_3+k_4) +A^{(3)}_{1,2,3,4}b_{k_2}^*b_{k_3}^*b_{k_4}^*\delta(k_1+k_2+k_3+k_4)\bigg].\end{align}

The goal here is that when (2.7) is substituted into (2.4), the resulting equation becomes

(2.8)\begin{equation}i{\dot b}_{k_1}=\omega_{k_1}b_{k_1}+\beta\sum_{k_1+k_2-k_3-k_4=0}T_{1,-2,3,4}b_{k_2}^\ast b_{k_3}b_{k_4}+O\left(\beta^2\right).\end{equation}

In order for (2.8) to be achieved, one needs to choose the coefficients of $A^{(i)}_{1,2,3,4}$ accordingly as

(2.9)\begin{equation} A^{(1)}_{1,2,3,4}=-\frac{T_{1,2,3,4}}{\omega_{k_1}-\omega_{k_2}-\omega_{k_3}-\omega_{k_4}}, \end{equation}

(2.10)\begin{equation} A^{(2)}_{1,2,3,4}=-\frac{T_{1,-2,3,4}}{\omega_{k_1}-\omega_{k_2}+\omega_{k_3}+\omega_{k_4}}, \end{equation}

(2.11)\begin{equation} A^{(3)}_{1,2,3,4}=-\frac{T_{1,2,3,4}}{\omega_{k_1}+\omega_{k_2}+\omega_{k_3}+\omega_{k_4}}. \end{equation}

3. Validity condition of the normal form transformation

In order for the normal form transformation (2.7) to stay bounded and near-identity, we need

(3.1)\begin{equation} \frac{\beta}{3}\left|\sum_{k_2,k_3,k_4}A^{(1)}_{1,2,3,4}b_{k_2}b_{k_3}b_{k_4}\delta(k_1-k_2-k_3-k_4)\right| \ll \left|b_{k_1}\right| \end{equation}

and similar conditions for other sums involving $A^{(2)}_{1,2,3,4}$ and $A^{(3)}_{1,2,3,4}$. A numerical study for similar conditions has been presented in [Reference Comito, Lotriglia and Bustamante4]. In addition, such a validity condition has also been formulated in [Reference Zakharov30] for gravity waves in finite water depth. However, in further estimating the sums in the LHS of the condition, the random phases in $b_k$ are not considered in [Reference Zakharov30], leading to an overestimation.

In order to account for the random phases in $b_k$, we consider a field with $b_k=\sqrt{\phi_k}\eta_k(\varrho)$, where $\phi_k$’s are deterministic amplitude of $O(1)$ and $\eta_k(\varrho)$’s are i.i.d. zero-mean complex random variables with unit variance, either as standard complex Gaussian or uniformly distributed on the unit circle. Under this condition, the sum in (3.1) behaves like a random walk in a complex plane, for which cancellations due to walks in different directions must be considered. Intuitively, this sum should be evaluated as the expected value $\mathbb{E}[|\sum \boldsymbol \cdot|]$ or $\{\mathbb{E}[|\sum \boldsymbol \cdot|^2] \}^{1/2}$. Mathematically, this reasoning is guaranteed by the standard hyper-contractivity estimates as in [Reference Deng and Hani6, Lemma 3.1], which asserts that the probability of the sum

(3.2)\begin{equation} \mathbb{P}\left(\left|\sum_{k_1-k_2-k_3-k_4=0}A^{(1)}_{1,2,3,4}b_{k_2}b_{k_3}b_{k_4}\right|\geq AM^{\frac{1}{2}}\right) \leq Ce^{-cA^{\frac{2}{3}}}, \end{equation}

for some positive constants $A,C,c \gt 0$, with

(3.3)\begin{equation} M=\mathbb{E}\left|\sum_{k_1-k_2-k_3-k_4=0}A^{(1)}_{1,2,3,4}b_{k_2}b_{k_3}b_{k_4}\right|^2. \end{equation}

This hyper-contractivity lemma simply states that with dominant probability, the sum in (3.1) is estimated by $M^{1/2}$.

By further applying the Wick contraction theorem (see [Reference Nazarenko20, Theorem 5.1]), we obtain

(3.4)\begin{equation} M = \sum_{k_1-k_2-k_3-k_4=0}\left|A^{(1)}_{1,2,3,4}\right|^2\phi_2\phi_3\phi_4. \end{equation}

The above calculation relies on the statistical independence of the modes. While this is assumed for the initial data, the persistence of such decorrelation, to leading order, under the non-linear flow is a non-trivial property known as propagation of chaos. As shown by Deng and Hani [Reference Deng and Hani10], for the non-linear Schrödinger (NLS) equation, the Wick contraction rules can be thought of as being preserved to leading order, in the kinetic limit (in our setting $N\to \infty$, $\beta\to 0$). More precisely, there holds that any moments of the flow of the form

(3.5)\begin{equation} \mathbb E \; b_{k_1} b_{k_2}\ldots b_{k_q}\overline{b_{\ell_1}}\, \overline{b_{\ell_2}}\ldots \overline{b_{\ell_r}} \to 0\end{equation}

in the kinetic limit, unless $\{k_1, \ldots, k_q\}=\{\ell_1, \ldots, \ell_r\}$. While the proof in [Reference Deng and Hani10] is done for the NLS equation, it is shown to be a consequence of the same diagrammatic analysis that is used to derive the WKE. It is shown that interactions that do not respect (3.5) are vanishing in the limit, and the leading order interactions, which appear in the kinetic limit, respect (3.5). This is yet another indication of two important facts: (a) Propagation of chaos is a rather deep property that cannot be simply explained by the linear theory; it is a deep non-linear effect; (b) propagation of chaos is a fundamental constituent of the kinetic theory, and that one should not expect the latter to hold without the former. As such, propagation of chaos will hold in much wider generality, namely for any system where the kinetic theory holds, in particular for the $\beta$-FPUT system considered here, for which the WKE has been rigorously derived [Reference Vassilev and Wu26, Reference Wu27]. This is the foundation for the scaling $M^{1/2}$ used here.

Finally, we caution that while the propagation of chaos property is sometimes referred to as preservation of the random phase property, phase randomness is not really preserved at the level of the dynamics of the microscopic system but is only valid as a limiting statement (i.e., the left hand side of (3.5) is not zero at positive times, but its limit is zero).

In the Appendix A of the paper, we prove the following theorem.

Theorem 3.1. Let $A^{(i)}_{1,2,3,4}$, $i=1,2,3$, be defined as in (2.9),(2.10) and (2.11), then for fixed $k_1 \in \mathbb{Z} \cap (0,N)$, we have:

(3.6)\begin{aligned}&\sum_{\substack{ k_1-k_2-k_3-k_4\equiv 0 \pmod N\\ k_2,k_3,k_4\in\mathbb{Z}\cap(0,N) }}\left|A_{1,2,3,4}^{(1)}\right|^2 +\sum_{\substack{ k_1-k_2+k_3+k_4\equiv 0 \pmod N\\ k_2,k_3,k_4\in\mathbb{Z}\cap(0,N) }}\left|A_{1,2,3,4}^{(2)}\right|^2 \\[6pt] &\quad +\sum_{\substack{ k_1+k_2+k_3+k_4\equiv 0 \pmod N\\ k_2,k_3,k_4\in\mathbb{Z}\cap(0,N)}} \left|A_{1,2,3,4}^{(3)}\right|^2 \lesssim N^2\log N. \end{aligned}

Combining (3.1) (and the counterparts for other sums), (3.4) and Theorem 3.1, we conclude that the validity condition of the normal form transformation is

(3.7)\begin{equation} \beta \ll \frac{1}{N^{1+\epsilon}}, \end{equation}

where $\epsilon$ is an arbitrarily small constant, resulting from the $\log N$ factor in (3.6) for large $N$. We remark that the sharp condition following from (3.1), (3.4) and Theorem 3.1 is $\beta \ll 1/(N\sqrt{\log N})$, with the correspondingly sharp controlling parameter $\beta N \sqrt{\log N}$. We adopt the more concise form (3.7) for two reasons. First, since $\log N$ grows more slowly than any positive power of $N$, the form (3.7) with $\epsilon$ arbitrarily small absorbs the logarithmic factor. Second, the logarithmic correction is too weak to be discerned over any practical range of system sizes, so that the numerical results cannot distinguish $\beta N \sqrt{\log N}$ from $\beta N$ as the controlling parameter, and $\beta N$ serves as the practical choice, consistent with the numerical observations in Section 4.

We note that to obtain (3.7), we have considered that $\phi_k \sim O(1)$ for all $k$. This means physically that the spectrum is broadband and does not deviate too much from the thermal-equilibrium state of equipartition of wave action. Of course, this is mathematically done for the convenience of cancellation of $\phi_i$’s in the condition. We will show numerically in next section that the condition derived under such a consideration does practically apply to spectra of different forms.

We remark that condition (3.7) admits an equivalent form that does not single out $\beta$. By the invariance of (2.1) under $q \to \lambda q$, $\beta \to \beta/\lambda^{2}$ (with $\lambda \gt 0$ an arbitrary amplitude scaling factor), the dynamics depends on $\beta$ and the amplitudes only through the products $\beta\phi_k$. The complete form of the validity condition thus reads $\beta\phi \ll 1/N^{1+\epsilon}$, with $\phi$ the common order of magnitude of the $\phi_k$’s, or equivalently $\beta e \ll 1/N^{1+\epsilon}$ (i.e., $\beta E \ll 1/N^{\epsilon}$) in terms of the energy per mode $e = E/N \propto \phi$, with $E = \sum_k \omega_k \phi_k$ the linear energy of the system. Condition (3.7) corresponds to the normalization $\phi_k \sim O(1)$, while setting $\beta = 1$ turns it into a condition on the energy $e \ll 1/N^{1+\epsilon}$. The condition is compatible with the kinetic limit of WT, i.e., the joint limit $N \to \infty$, $\beta \to 0$ under a scaling law between $\beta$ and $N$ such that (3.7) remains satisfied along the limit. In particular, the scaling laws under which the WKE is rigorously derived for the FPUT system [Reference Vassilev and Wu26, Reference Wu27] lie within this admissible regime (with the non-linearity parameter normalized consistently; see also [Reference Deng and Hani5–Reference Deng and Hani9] for NLS), so that the normal form transformation remains valid along the kinetic limit.

4. Numerical verification

The condition (3.7) suggests that the relative importance of the three non-resonant terms in (2.4) to the dynamics is governed by the parameter $\beta N$ ( $\epsilon$ omitted hereafter), which we aim to verify numerically in this section. We emphasize that the purpose of the numerical study is not a direct verification of condition (3.1), which is established rigorously under the random-phase assumption through Theorem 3.1. Instead, we verify a central physical consequence of the theory: that the relative contribution of the non-resonant interactions, measured in terms of the energy content (Hamiltonian) of the system, is governed by the single parameter $\beta N$, irrespective of the values of $\beta$ and $N$ separately. We first note that the Hamiltonian for the original $\beta$-FPUT system is given by

(4.1)\begin{equation} \begin{split} &\frac{H}{N}=\sum_{k=1}^{N-1}\omega_k|a_k|^2+ {\frac{\beta}{3}}\sum_{k_1,k_2,k_3,k_4} \bigg[ T_{1,2,3,4}(a_1^*a_2 a_3 a_4+c.c.) \delta_{1-2-3-4} \\ &\qquad\qquad\qquad\qquad\qquad+ \frac{3}{2}T_{1,-2,3,4}a_1^*a_2^* a_3 a_4 \delta_{1+2-3-4}+\frac{1}{4}T_{1,2,3,4}(a_{1}^*a_{2}^*a_{3}^*a_{4}^*+c.c. )\delta_{1+2+3+4} \bigg], \end{split} \end{equation}

where $c.c.$ denotes the complex conjugate of the previous term. The relative dynamical contribution from the non-resonant terms can be characterized by the ratio of non-resonant to resonant terms in (4.1). More specifically, we first define

(4.2)\begin{equation} S_{1}=\sum_{k_1,k_2,k_3,k_4} T_{1,2,3,4}(a_1^*a_2 a_3 a_4+c.c.) \delta_{1-2-3-4}, \end{equation}

(4.3)\begin{equation} \!\!\!\!\!\!\!\!\!\!\!\!S_{2}=\sum_{k_1,k_2,k_3,k_4}\frac{3}{2}T_{1,-2,3,4}a_1^*a_2^* a_3 a_4 \delta_{1+2-3-4}, \end{equation}

(4.4)\begin{equation} S_{3}=\sum_{k_1,k_2,k_3,k_4}\frac{1}{4}T_{1,2,3,4}(a_{1}^*a_{2}^*a_{3}^*a_{4}^*+c.c. )\delta_{1+2+3+4}, \end{equation}

and then the ratio

(4.5)\begin{equation} r = \left|\frac{\overline{S_1}+\overline{S_3}}{\overline{S_2}}\right|, \end{equation}

where $\overline{\boldsymbol \cdot}$ represents average of the quantity (computed from time and ensemble average specified later). Our goal is to numerically compute $r$ as a function of $\beta$ and $N$ and show that $r$ depends on $\beta N$ irrespective of each factor, as theory suggests. We note that the ratio $r$ compares the energy content of the non-resonant and resonant parts of the quartic Hamiltonian, rather than their contributions to the energy transfer. In particular, the diagonal terms in (4.3) (i.e., the trivial resonances with $\{k_1, k_2\} = \{k_3, k_4\}$) produce the non-linear frequency correction and do not contribute to the spectral energy transfer described by the WKE. They are nevertheless part of the resonant Hamiltonian retained by the normal form transformation, entering the dynamics through the frequency renormalization, and are therefore included in $S_2$, consistent with our purpose of quantifying the fraction of the quartic Hamiltonian carried by the non-resonant terms relative to the resonant terms retained after the transformation.

For the purpose mentioned above, we perform numerical simulations of (2.1) using the sixth order symplectic integrator method [Reference Yoshida28]. We choose two sets of different initial conditions, one at thermal equilibrium (of energy equipartition):

(4.6)\begin{equation} a_k{(0)}=\sqrt{1/\omega_k}\cdot\exp(i\varphi_k), \end{equation}

and one out of equilibrium:

(4.7)\begin{equation} a_k{(0)}=\sqrt{(1+\omega_k^2)/\omega_k}\cdot\exp(i\varphi_k), \end{equation}

with phases $\varphi_k$ uniformly distributed from $[0,2\pi)$ for all $k\in \mathbb{Z}\cap(0,N)$. We can then recover $q_j(0),p_j(0)$ from (2.2) and (2.3), as the initial condition for (2.1). We set $m = 1, \kappa = 1$, so that the fundamental period $T_f=\frac{2\pi}{\omega(k=1)}\approx N$. We set the time step $h=0.01$ and simulate the system up to $t_{\text{max}}=10T_f$. The average in (4.5) is computed from time average over $[5T_f, 10T_f]$, then ensemble average over 5 ensembles.

Figure 1.

Ratio $r$ as a function of $\beta N$ for $N = 200$, $500$, $800$ and $1000$, with initial conditions (a) at thermal equilibrium with (4.6) and (b) out of thermal equilibrium with (4.7).

1 Long description

Two scatter plots labeled (a) and (b). The horizontal axis in both plots is labeled beta N, ranging from 0.01 to 1. The vertical axis is labeled r, ranging from 0 to 0.90. In plot (a), titled Initial Data At Equilibrium, the data points are plotted for different system sizes: N equals 200, N equals 500, N equals 800 and N equals 1000. The legend indicates different markers for each system size: plus for N equals 200, cross for N equals 500, diamond for N equals 800 and circle for N equals 1000. The data points show a trend where r increases with beta N, with values approximately at (0.01, 0.04), (0.05, 0.16), (0.1, 0.27), (0.5, 0.64) and (1, 0.80). In plot (b), titled Initial Data Out of Equilibrium, the same system sizes and markers are used. The trend is similar to plot (a), with r increasing as beta N increases. The data points are approximately at (0.01, 0.06), (0.05, 0.22), (0.1, 0.35), (0.5, 0.74) and (1, 0.87). Both plots demonstrate that the values of r collapse for the same value of beta N irrespective of the system size.

Figure 1 plots the computed values of $r$ from the two initial conditions (4.6) and (4.7) for different values of $\beta$ and $N$. In particular, we choose $N = 200$, $500$, $800$ and $1000$, and for each $N$ corresponding values of $\beta$ to make $\beta N=1$, $0.5$, $0.1$, $0.05$ and $0.01$. We see from Figure 1 that all results of $r$ collapse for the same value of $\beta N$ irrespective of the value of $N$, which is true for both initial conditions. This means that the relative contributions from non-resonant terms are indeed controlled by the parameter $\beta N$ as indicated by the theory. In addition, we see that for large values of $\beta N$, say $\beta N=1$, the contribution from non-resonant terms can rise up to $O(80\%)$ for both initial conditions. Clearly for these cases, ignoring the non-resonant terms (i.e., applying the normal form transformation inconsistently) leads to excessive errors to the dynamics of the system.