Anomalous correlators,"ghost"waves and nonlinear standing waves in the $\beta$-FPUT system

We show that Hamiltonian nonlinear dispersive wave systems with cubic nonlinearity and random initial data develop, during their evolution, anomalous correlators. These are responsible for the appearance of"ghost'' excitations, i.e. those characterized by negative frequencies, in addition to the positive ones predicted by the linear dispersion relation. We explain theoretically the existence of anomalous correlators by nonresonant wave-wave interactions, extending {some} concepts widely used in the Wave Turbulence Theory. Namely, by generalizing the classical Wick's decomposition to include the second-order anomalous correlator, we show that the latter is responsible for the presence of ``ghost'' excitations. We test our theory on the celebrated $\beta$-Fermi-Pasta-Ulam-Tsingou chain and show that numerically measured values of the anomalous correlators agrees, in the weakly nonlinear regime, with our analytical predictions. We also predict that similar phenomenon will occur in other nonlinear systems dominated by nonlinear interactions, including surface gravity waves. Our results paves the road to study phase correlations in the Fourier space.


I. INTRODUCTION
Wave Turbulence theory has led to successful predictions on the wave spectrum in many fields of physics [1,2]. In this framework the system is represented as a superposition of a large number of weakly interacting waves with the complex normal variables a k = a(k, t). In its essence, the classical Wave Turbulence theory is a perturbation expansion in the amplitude a k of the nonlinearity, yielding, at the leading order, to a system of quasilinear waves whose amplitudes are slowly modulated by resonant nonlinear interactions [3][4][5][6][7]. This modulation leads to a redistribution of the spectral energy density among length-scales, and is described by a wave kinetic equation. One way to derive the wave kinetic equation is to use the random phase and amplitude approach developed in [2,8,9]. The initial state of the system can be always prepared so that the assumption of random phases and amplitudes is true. Whether the phases remain random in the evolution of the system has been an issue of intense discussions. In Wave Turbulence theory, the standard object to look at is the second-order correlator, a k (t)a * l (t) , where . . . is an average over an ensemble of initial conditions with different random phases and amplitudes. As will be clear later on, under the homogeneity assumption, the second order correlator is related to the wave action spectral density function, i.e. the wave spectrum, n k = n(k, t). However, one should note that the complex normal variable, as defined in the the Wave Turbulence theory, is a complex function also in physical space. Therefore, the second-order statistics are not fully determined by the above correlator. The so called "anomalous correlator", a k (t)a l (t) , see [10,11], needs also to be computed. Under the hypothesis of homogeneity, will be related to the anomalous spectrum, m k = m(k, t), to be defined in the next Section. Indeed, if phases are totally random, this quantity would be zero. We show that, in the nonlinear evolution of the system, this is not the case. Far from it, this quantity is strongly nonzero and, in the limit of weak nonlinearity, we predict analytically and verify numerically its value.
Our ideas are based on the extension of the Wave Turbulence Theory which has been successfull in the understanding of the spectral energy transfer in complex wave systems such as the ocean [12], opticas [13] and Bose-Einstein condensates [14], one dimensional chains [15], and magnets [16]. One of the main tool used to derive the theory is the Wick's contraction rule that allows one to split higher order correlators as sum of products of second order ones, plus cumulants. We show that the rule used is incomplete and needs an extra terms, namely the products of anomalous correlators, to the presence of "ghost waves" observed in numerical simulations.
These ideas are tested on a simple, but non trivial, system, i.e. the celebrated β-Fermi-Pasta-Ulam-Tsingou (FPUT) chain. The chain model was introduced in the fifties to study the thermal equipartition in crystals [17]; it consists of N identical masses each one connected by a nonlinear spring; the elastic force can be expressed as a power series in the spring deformation. Fermi, Pasta, Ulam and Tsingou integrated numerically the equations of motion and conjectured that, after many iterations, the system would exhibit a thermalization, i.e. a state in which the influence of the initial modes disappears and the system becomes random, with all modes excited equally (equipartition of energy) on average. Successful predictions on the time scale of equipartition have been recently obtained in [18][19][20] using the Wave Turbulence approach. In this paper we perform extensive numerical simulations with initial random data and look all at the possible excitations, once a thermalized state has been reached. This is all done by analysing the spatial -temporal (k − Ω) spectrum, i.e. the square of the spacetime Fourier transform of the wave amplitudes. Analyses of the effective dispersion relation in the nonlinear system is a well known and widely used theoretical and numerical tool, see for example [21].
We give numerical evidence that in addition to the "normal" waves with frequency ω predicted by the linear dispersion relation for wave number k, there are the "ghost" excitations with the negative frequencies.
Our theoretical analysis reveals that the origin of those "ghost" excitations resides on the nonzero values of the second-order anomalous correlator.

II. THE MODEL
The theory that we develop hereafter applies to any system with cubic nonlinearity. Examples of such systems among others, include surface gravity waves on top of a deep fluid [1], Nonlinear Schrodinger equation, β-Fermi-Pasta-Ulam-Tsingou chain and Majda-McLaughlin-Tabak system [22]. In normal variables a k the Hamiltonian of these systems assumes the canonical form: where ω k = ω(k) are the frequencies associated to the wave numbers via the dispersion relation, T 1234 are coefficients that depend on the problem considered and satisfy specific symmetries for the system to be Hamiltonian, c.c. implies complex conjugation, a j = a(k j , t) are the complex normal variables, δ lm ij = δ(k i + k j − k l − k m ) is the Kronecker Delta. We assume that the only resonant interactions possible are the ones for which the following two relations are satisfied for a set of wave numbers With the objective of presenting some comparison with numerical simulations, out of many physical systems described by the above Hamiltonian, we select a simple one dimensional system, i.e. the Fermi-Pasta-Ulam-Tsingou chain. Therefore, our problem consists of a system of N identical particles connected locally to each other by a nonlinear spring. In the physical space the displacements with respect to the equilibrium position q j (t) and their momenta p j (t), the Hamiltonian takes the following form: with β is the nonlinear spring coefficient (without loss of generality, we have set the masses and the linear spring constant equal to 1). The Newton's law in physical space is given by: We assume periodic boundary condition; our approach is developed in Fourier space and the following definitions of the direct and inverse Discrete Fourier Transforms are adopted: where k are discrete wave numbers and Q k are the Fourier amplitudes. The displacement q j and momentum p j of the j particle are linked by canonically conjugated Hamilton equationsṗ We then perform the Fourier transformation to Q k , P k , and then additional canonical transformation to complex amplitude a k given by where ω k = 2| sin(πk/N )| > 0 and Q k and P k are the Fourier amplitudes of q j and p j , respectively. In terms of a k the equation of motion reads: where all wave numbers k 2 , k 3 and k 4 are summed from 0 to N − 1 and δ cd.. ab.
is the generalized Kronecker Delta that accounts for a periodic Fourier space, i.e. its value is one when the argument is equal to 0 (mod N ). Furthermore, The main statistical object discussed in this paper is the wave number-frequency (k − Ω) spectrum. Starting from the complex amplitude a(k, t) we take the Fourier transform in time so that we get a(k, Ω); Under the hypothesis of homogeneous and stationary conditions, the second-order (k − Ω) correlator takes the following form where .... implies averages over initial conditions with different random phases. N (k, Ω) is the (k −Ω) spectrum defined as follows: The linear (k −Ω) spectrum -Before diving into the nonlinear dynamics, we discuss the predictions in the linear regime. Therefore, we start by neglecting the nonlinearities in equation (8) and find the solution in the form where t 0 is a time at which the solution is known or an initial condition. We then take Fourier transform in time After multiplication by its complex conjugate and after taking averages, we get: where n (a) (k, t 0 ) is the standard wave spectrum at time t 0 related to the second-order correlator as and defined via the autocorrelation function as In the linear regime n (a) (k i , t 0 ) does not evolve in time. Equation (14) implies that in the linear case the (k−Ω) spectrum is different from zero only for those values of Ω and k for which the dispersion relation is satisfied. Note that in this formulation ω k is defined as a positive quantity; therefore, only the positive branch of the dispersion relation curve appears in the linear regime. We now test the predictions from equation (14) both in the linear regime and observe what happens to it in the nonlinear regime. We perform numerical simulations of the equations (5) using a symplectic algorithm, see [23]. We use 32 particles in the simulations; such choice is completely uninfluential for the results presented below. In the linear regime, we just prescribe a thermalized spectrum with some initial random phases of the wave amplitudes a k and evolve the system in time up to a desired final time; a Fourier Transform in time is then taken to build the (k − Ω) spectrum. In the nonlinear regime we perform long simulations up to a thermalized spectrum. For a given nonlinearity, about 1000 realizations characterized by different random phases are made and ensemble averages are considered to compute the (k − Ω) spectrum. All simulations have the same initial linear energy and, from an operative point of view, the only difference between them is the value of β. To charactrerize the strength of the nonlinearity, we use the following ratio between nonlinear and linear Hamiltonians at the beginning of each simulation: Results are shown in Figure 1, where, for different values of the nonlinear parameter , the spectrum N (a) (k, Ω) is plotted using a coloured logarithmic scale. We first focus our attention on the linear regime, = 0: results are shown in Figure 1(a). As well predicted by the theory, the plot shows dots in the positive frequency plane, where the frequencies Ω and wave numbers k satisfy the linear dispersion curve ω k . Increasing the nonlinearity, Figures 1 -(b,c,d), two well known effects appears: the first one is a shift of the frequencies, due to nonlinearity (this is more evident in Figures 1 -(c,d) where the frequency scale in the vertical axes has been changed). The second one is the broadening of the frequencies; this is related to the fact that the amplitude for each wave number is not constant in time; therefore, the amplitudedependent frequencies are not constant in time and they oscillate around a mean value with some fluctuations. Those results are well understood, at least in the weakly nonlinear regime, and can be predicted using Wave Turbulence tools, see [2,20]. Besides these two effects, starting from Figure 1-(b), the presence of a lower branch, whose intensity is much less than the upper one, starts to be visible. The lower curve becomes more important and, when the nonlinearity is of order one, is of the same order of magnitude of the upper one. The total number of waves in the simulation, N tot , is given by the integral over Ω and the sum over all k of the function N (a) (k, Ω). In the weakly nonlinear regime, N tot is an adiabatic invariant of the equation of motion (5); the plot highlights the existence of waves with negative frequencies, which will be named "ghost" exitations. One of the scopes of the present paper is the understanding of the origin of such waves. Before entering into the discussion, we show in Figure II B the number of "ghost" excitations, N ghost , i.e. N (a) (k, Ω) integrated over negative frequencies and summed over all wave numbers, divided by the total number of waves, N tot . As can be seen from the plot, there is a monotonic growth of the "ghost" waves that, for very large nonlinearity, can reach values up to 25% of the total number.

III. ANOMALOUS CORRELATORS
To explain the presence of "ghost" excitations, we introduce the so called second-order anomalous correlator [10,11,24]: with the anomalous spectrum defined as m (a) Similarly, we also introduce the second-order (k − Ω) anomalous correlator:  To verify numerically that the anomalous correlator is indeed nonzero, we measure numerically the real part of the second-order correlator a ki (t)a kj (t) as a function of k 1 and k 2 . Results are plotted in Figure III  (a) = 0.0089 and (b) = 1.12. In both cases, a diagonal contribution is visible, pointing out the existence of anomalous correlators in the β-FPUT model. Generalization of the Wick's decomposition -It is now straightforward to extend the Wick's decomposition by taking into account the anomalous correlators as follows [16]: a * k a * l a p a n = n k n l (δ k p δ l n + δ k n δ l p ) + m * k m p δ kl δ pn , a * k a l a p a n = n k m p (δ l k δ pn + δ n k δ lp ) + n k δ p k m l δ nl , a k a l a p a n = m k m l (δ kp δ ln + δ kl δ pn + δ kl δ pn ). (22) The above relations will be fundamental for making a natural closure of the moments when calculating analytically the (k − Ω) spectrum.
We show in Figure V the real part of the fourth-order correlator a k1 a k2 a * k3 a * k1+k2−k3 with k 3 = 20, computed from numerical simulations for (a) = 0.0089 and (b) = 1.12. The diagonal lines in both figures, highlighting the contribution from the second-order anomalous correlator, are noticeable. The vertical and horizontal lines correspond to the trivial resonances in which two wave numbers are equal (mod N ).

A. Theoretical prediction for the anomalous correlator in the weakly nonlinear regime
A key step for the development of a theory for the anomalous correlator is the change of variable (near identity transformation) which allows one to remove bound modes, i.e. those modes that are phase locked to free modes and do not obey the linear dispersion relation. The procedure is well known in Hamiltonian mechanics and well documented for example in [1]. The transforma-tion from variable a k (t) to b k (t) has the following form: The coefficients B (i) 1234 are selected in such a way to remove non resonant terms in the original Hamiltonian. The transformation is asymptotic in the sense that the small amplitude approximation is made and the terms in the sums on the right hand side are much smaller than the leading order term b 1 . The evolution equation for variable b k (t) contains resonant interactions and take the following standard form: where higher order terms arising from the transformation have been neglected. Starting from the transformation in (23) and the generalized Wick's decomposition in (22), we can now build the anomalous spectrum: where higher order terms in m k have been neglected. The next step consists in building the evolution equation for m  k (t) appears as a deterministic dispersive non homogeneous wave evolution equation The leading order solution is given by k (t 0 )e −i2ω k t + higher order terms.
Putting things together, and assuming that the spectrum n k is in stationary conditions, we get m (a) Note that we have used the fact that at the leading order n k (t 0 ). Averaging over time, in the limit of large times, assuming that the wave spectrum n where ... t implies averaging over time. The prediction is shown in figures 5 against numerical simulations: the results are in good agreement in the weakly nonlinear regime.

IV. THEORETICAL PREDICTION FOR "GHOST" EXCITATIONS
We have now developed all the tools for predicting analytically the (k −Ω) spectrum as defined in the equation (10). We start by writing the transformation in equation (23) in the frequency domain as follows: and build the correlator a(k i , Ω l )a(k j , Ω m ) * . Using the generalized Wick's decomposition and the hypothesis of stationarity and homogeneity, we get at leading order: where we have used the fact that at the leading order m (b) (k, t) m (a) (k, t) and n (b) (k, t 0 ) n (a) (k, t 0 ). The   (31) predict the presence of the upper and lower branch in the (k − Ω) plane. The presence of "ghost" excitations is clearly related to the second-order anomalous correlator. We can now predict the percentage of "ghost" excitations as In figure (6) we plot the ratio (32) for several values of . We find that for small nonlinearity the results agree.

V. NONLINEAR STANDING WAVES
The development of a regime characterized by an anomalous spectrum m (a) k (t) t corresponds to a tendency for the system to develop standing waves in the original displacement q j (t). Indeed, the existence of an an anomalous spectrum implies a correlation between positive and negative wave numbers. We illustrate such behavior by initializing the system in Fourier space with initial conditions with all particles at rest and displaced from the equilibrium as a single sine wave. In figure (7) we plot the displacement for all masses as a function of time as system reaches thermal equilibrium. The existence of several regions of standing wave behavior is clearly visible.

VI. CONCLUSION
In this paper we have developed a theory for weakly nonlinear dispersive waves that accounts for presence of an anomalous correlator. The framework in which the theory has been developed is the Wave Turbulence one. In such theory one usually is interested in the second order correlator a ki a * kj which is strictly related to the wave action spectrum. However, what is clear from numerical simulations of the β-FPUT system is that also the correlator a ki a kj can assume values that are different from zero. This finding has consequences on the standard Wave Turbulence theory that is based on the Wick's selection rule, i.e. the splitting of higher order correlators as a sum of products of second order correlators. Following [1], we have generalized the Wick's rule by including the anomalous correlators. One of the most striking manifestation of those correlators is the appearance of "ghost excitations", i.e. those characterized by a negative frequencies. A formula for the content of energy of such excitations as a function of the wave spectrum is obtained and compared favourably, in the weakly nonlinear, regime with numerical simulations.
Our approach paves a new road to investigate dispersive nonlinear systems by taking into account not only amplitudes of the waves, as in traditional wave turbulence, but also the phases of the waves. We conjecture that the anomalous correlators play an important role in the theory of extreme events such as for example rogue waves, [25] and in addressing the statistical properties of integrable systems such as the Nonlinear Schrodinger equation. Moreover, being the Hamiltonian here considered of the same family as the one for surface gravity waves (after removing by a canonical transformation nonreaonant three wave interactions), we predict that also the anomalous correlators may play an important role in the understanding of statistical properties of ocean waves.