Numerical study of the $\beta$-FPU chain

Since its introduction in the early 1950s, the study of the FPU lattice [13] has led to many great discoveries in mathematics and physics, such as soliton theory [3]. Being non-integrable, the FPU system also became intertwined with the celebrated Kolmogorov-Arnold-Moser theorem [11]. Here, we extend our results of the thermalized $\beta$-FPU chain, which were briefly reported in [12].

The Hamiltonian of the $\beta$-FPU chain is of the form

$$H=\sum_{j=1}^N\frac{1}{2}p_j^2+\frac{1}{2}(q_j-q_{j+1})^2+\frac{\beta}{4}(q_j-q_{j+1})^4,$$ (29)

where $\beta$ is a parameter that characterizes the strength of nonlinearity.

The canonical equations of motion of the $\beta$-FPU chain are

$$\dot{q}_j = \frac{\partial H}{\partial p_j}=p_j,$$ (30)

$$\dot{p}_j = -\frac{\partial H}{\partial q_j}=(q_{j-1}-2q_j+q_{j+1})$$

$$+ \beta\big[(q_{j+1}-q_j)^3-(q_j-q_{j-1})^3\big].$$

To investigate the dynamical manifestation of the renormalized dispersion $\tilde{\omega}_k$ of $\tilde{a}_k$, we numerically integrate Eq. (30). Since we study the thermal equilibrium state [14,15,16,17] of the $\beta$-FPU chain, we use random initial conditions, i.e., $p_j$ and $q_j$ are selected at random from the uniform distribution in the intervals $(-p_{\rm {max}},p_{\rm {max}})$ and $(-q_{\rm {max}},q_{\rm {max}})$, respectively, with the two constraints that (i) the total momentum of the system is zero and (ii) the total energy of the system $E$ is set to be a specified constant. We have verified that the results discussed in the paper do not depend on details of the initial data. Note that the behavior of $\beta$-FPU for fixed $N$ is fully characterized by only one parameter $\beta E$ [18]. We use the sixth order symplectic Yoshida algorithm [19] with the time step $dt=0.01$, which ensures the conservation of the total system energy up to the ninth significant digit for a runtime $\tau=10^6$ time units. In order to confirm that the system has reached the thermal equilibrium state [20], the value of the energy localization [21] was monitored via $L(t)\equiv{N\sum_{j=1}^{N}G_j^2}/{(\sum_{j=1}^{N}G_j)^2}$, where $G_j$ is the energy of the $j$-th particle defined as

$$G_j = \frac{1}{2}p_j^2+\frac{1}{4}\big[(q_j-q_{j+1})^2+(q_{j-1}-q_j)^2\big]$$

$$+ \frac{\beta}{8}\big[(q_j-q_{j+1})^4+(q_{j-1}-q_j)^4\big].$$ (31)

If the energy of the system is concentrated around one site, then $L(t)=O(N)$. Whereas, if the energy is uniformly distributed along the chain, then $L(t)=O(1)$. In our simulations, in thermal equilibrium states, $L(t)$ is fluctuating in the range of $1$-$3$. Since our simulation is of microcanonical ensemble, we have monitored various statistics of the system to verify that the thermal equilibrium state that is consistent with the Gibbs distribution (canonical ensemble) has been reached. Moreover, we verified that, for $N$ as small as 32 and up to as large as 1024, the equilibrium distribution in the thermalized state in our microcanonical ensemble simulation is consistent with the Gibbs measure. We compared the renormalization factor (25) by computing the values of $\langle K\rangle$ and $\langle U\rangle$ numerically and theoretically using the Gibbs measure and found the discrepancy of $\eta$ to be within $0.1\%$ for $\beta=1$ and the energy density $E/N=0.5$ for $N$ from 32 to 1024.

We now address numerically how the renormalized linear dispersion $\tilde{\omega}_k$ manifests itself in the dynamics of the $\beta$-FPU system. We compute the spatiotemporal spectrum $\vert\hat{a}_k(\omega )\vert^2$, where $\hat{a}_k(\omega )$ is the Fourier transform of $\tilde{a}_k(t)$. (Note that, for simplicity of notation, we drop a tilde in $\hat{a}_k$.)

Figure: The spatiotemporal spectrum $\vert\hat{a}_k(\omega )\vert^2$ in thermal equilibrium. The chain was modeled for $N=256$, $\beta =0.5$, and $E=100$. [ $\max\{-8,\ln{\vert\hat{a}_k(\omega )\vert^2}\}$, with corresponding gray scale, is plotted for a clear presentation]. The solid curve corresponds to the usual linear dispersion $\omega _k=2\sin(\pi k/N)$. The dashed curve shows the locations of the actual frequency peaks of $\vert\hat{a}_k(\omega )\vert^2$.

Figure 1 displays the spatiotemporal spectrum of $\tilde{a}_k$, obtained from the simulation of the $\beta$-FPU chain for $N=256$, $\beta =0.5$, and $E=100$. In order to measure the value of $\eta$ from the spatiotemporal spectrum, we use the following procedure. For the fixed wave number $k$, the corresponding renormalization factor $\eta (k)$ is determined by the location of the center of the frequency spectrum $\vert\hat{a}_k(\omega )\vert^2$, i.e.,

$$\eta(k)=\frac{\omega _{c}(k)}{\omega _k},~~ ~~ \omega _{c}(k)=\frac{{\int} \omega \vert\hat{a}_k(\omega )\vert^2~d\omega }{\int\vert\hat{a}_k(\omega )\vert^2~d\omega }.$$

The renormalization factor $\eta (k)$ of each wave mode $k$ is shown in Fig. 2 (inset). The numerical approximation $\bar{\eta}$ to the value of $\eta$ is obtained by averaging all $\eta (k)$, i.e.,

$$\bar{\eta}=\frac{1}{N-1}\sum_{k=1}^{N-1}\eta(k).$$

The renormalization factor for the case shown in Fig. 1 is measured to be $\bar{\eta}\approx1.1824$. It can be clearly seen in Fig. 2 (inset) that $\eta (k)$ is nearly independent of $k$ and its variations around $\bar{\eta}$ are less than $0.3\%$. We also compare the renormalization factor $\eta$ obtained from Eq. (25) (solid line in Fig. 2 (inset)) with its numerically computed approximation $\bar{\eta}$ (dashed line in Fig. 2 (inset)). Equation (25) gives the value $\eta \approx 1.1812$ and the difference between $\eta$ and $\bar{\eta}$ is less then $0.1\%$, which can be attributed to the statistical errors in the numerical measurement.

Figure 2: The renormalization factor as a function of the nonlinearity strength $\beta$. The analytical prediction [Eq. (25)] is depicted with a solid line and the numerical measurement is shown with circles. The chain was modeled for $N=256$, and $E=100$. Inset: Independence of $k$ of the renormalization factor $\eta (k)$. The circles correspond to $\eta (k)$ obtained from the spatiotemporal spectrum shown in Fig. 1 [only even values of $k$ are shown for clarity of presentation]. The dashed line corresponds to the mean value $\bar{\eta}$. For $\beta =0.5$, the mean value of the renormalization factor is found to be $\bar{\eta}\approx1.1824$. The variations of $\eta _k$ around $\bar{\eta}$ are less then $0.3\%$. [Note the scale of the ordinate.] The solid line corresponds to the renormalization factor $\eta$ obtained from Eq. (25). For the given parameters $\eta \approx 1.1812$.

In Fig. 2, we plot the value of $\eta$ as a function of $\beta$ for the system with $N=256$ particles and the total energy $E=100$. The solid curve was obtained using Eq. (25) while the circles correspond to the value of $\eta$ determined via the numerical spectrum $\vert\hat{a}_k(\omega )\vert^2$ as discussed above. It can be observed that there is excellent agreement between the theoretic prediction and numerically measured values for a wide range of the nonlinearity strength $\beta$.

In the following Sections, we will discuss how the renormalization of the linear dispersion of the $\beta$-FPU chain in thermal equilibrium can be explained from the wave resonance point of view. In order to give a wave description of the $\beta$-FPU chain, we rewrite the Hamiltonian (29) in terms of the renormalized variables $\tilde{a}_k$ [Eq. (19)] with $\tilde{\omega}_k=\eta\omega _k$,

$$H = \sum_{k=1}^{N-1}\frac{\omega _k}{2}\left(\eta+\frac{1}{\eta}\right)\vert\tilde{a}_k\vert^2$$

$$+ \frac{\omega _k}{4}\left(\eta-\frac{1}{\eta}\right)(\tilde{a}_k^*\tilde{a}_{N-k}^*+\tilde{a}_k\tilde{a}_{N-k})$$

$$+ \sum_{k,l,m,s=1}^{N-1}T^{kl}_{ms}\Big[\Delta _{ms}^{kl}\tilde{a}_k\tilde{a}_l\tilde{a}_m^*\tilde{a}_s^*$$

$$+ \left(\frac{2}{3}\Delta ^{klm}_{s}\tilde{a}_k\tilde{a}_l\tilde{a}_m\tilde{a}_s^*+c.c.\right)$$

$$+ \left(\frac{1}{6}\Delta ^{klms}_{0}\tilde{a}_k\tilde{a}_l\tilde{a}_m\tilde{a}_s+c.c.\right)\Big],$$ (32)

where c.c. stands for complex conjugate, and

$$T^{kl}_{ms}=\frac{3\beta}{8N\eta^2}\sqrt{\omega _k\omega _l\omega _m\omega _s}$$ (33)

is the interaction tensor coefficient. Note that, due to the discrete nature of the system of finite size, the wave space is periodic and, therefore, the ``momentum'' conservation is guaranteed by the following ``periodic'' Kronecker delta functions

$$\Delta ^{kl}_{ms} \equiv \delta_{ms}^{kl}-\delta^{klN}_{ms}-\delta^{kl}_{msN},$$ (34)

$$\Delta ^{klm}_{s} \equiv \delta^{klm}_{s}-\delta^{klm}_{sN}+\delta^{klm}_{sNN},$$ (35)

$$\Delta ^{klms}_{0} \equiv \delta^{klms}_{NN}-\delta^{klms}_{N}-\delta^{klms}_{NNN}.$$ (36)

Here, the Kronecker $\delta$-function is equal to 1, if the sum of all superscripts is equal to the sum of all subscripts, and 0, otherwise.