Statistical analysis of nonlinear wave interactions in simulated Langmuir turbulence data

We present a statistical analysis of strong turbulence of Langmuir and ion-sound waves resulting from beam-plasma interaction. The analysis is carried out on data sets produced by a numerical simulation of one-dimensional Zakharov’s equations. The nonlinear wave interactions are studied using two different approaches: high-order spectra and Volterra models. These methods were applied to identify two and three wave processes in the data, and the Volterra model was furthermore employed to evaluate the direction and magnitude of energy transfer between the wave modes in the case of Langmuir wave decay. We demonstrate that these methods allow one to determine the relative importance of strongly and weakly turbulent processes. The statistical validity of the results was thoroughly tested using surrogated data set analysis.


Introduction
The understanding of Langmuir turbulence in space plasmas is a longstanding research topic that has received much attention both from a theoretical and from a numerical point of view (see Robinson, 1997;Goldman et al., 1996, for a review).The classical scenario for the generation of weak plasma turbulence generally admits three stages: 1) energetic electron beams generate Langmuir waves by means of the beam-plasma instability, 2) these primary waves transfer their energy to lower frequency secondary waves by nonlinear wave-wave interactions, 3) finally, wave-particle interactions stabilize the non-maxwellian particle distribution.The turbulence that is generated that way has traditionally been divided into weak and strong turbulence.
This scenario is subject to controversy, and the study of its solutions has proved to be a difficult task.Not only does Correspondence to: J. Soucek (soucek@ufa.cas.cz)one need to make many assumptions (such as the random phase approximation for the weak turbulence case) but, more importantly, the comparison of the solutions against experimental or numerical data has often remained inconclusive.One of the reasons for this is that we still lack numerical techniques that can extract the key features of Langmuir turbulence, and thereby provide the opportunity to test theory against experiment.
The objective of our study is to address this problem of the characterization of statistical properties of Langmuir turbulence.We'll focus on one particular aspect only, which is the parametric decay of primary waves by means of nonlinear wave-wave interactions.Indeed, it is well known that the weak turbulence regime can be adequately described in terms of three-wave and four-wave interactions (Galeev and Sudan, 1989).The appropriate techniques for quantifying the dynamics and the statistical characteristics of such interactions are higher order spectra and Volterra models.Here we shall use simulation data to show how these techniques work, what are their validity limits, and how they can provide a new insight into the underlying physics.It must be stressed that these same techniques are applicable to many other processes, whenever nonlinear wave interactions occur in a conservative system (Kadomtsev, 1965).Examples of such wave-wave interactions can be found in the ionosphere, in ocean waves, in nonlinear optics, in chemical reactions, etc.

Simulation data
Our study is based on a code that simulates 1-D Langmuir turbulence in a non-isothermal plasma (T e T i ), which is excited by an electron beam with a frequency above ω pe .The properties of this system have been thoroughly investigated, both theoretically and numerically (Shapiro and Shevchenko, 1983;Al'terkop et al., 1976).
The code integrates the Zakharov equations (Zakharov, 1972) in one spatial dimension (2) These equations describe the time evolution of highfrequency Langmuir oscillations, together with the lowfrequency fluctuations of plasma density.The electric field E appearing in these equations is a complex quantity related to the real physical field Ẽ as Ẽ(x, t) = 1 2 E(x, t)e −iω pe t + E * (x, t)e iω pe t (3) and ρ(x, t) = δn/n 0 describes the electron density fluctuations.The operators ν and γ S specify the electron and ion damping rates.The equations are Fourier transformed in space coordinate (using periodic boundary conditions) and then integrated using a standard finite difference scheme.
The simulation starts from an initial weak random noise and the energy necessary to develop the turbulence is delivered to the system by the oscillating electron beam.This is in practice achieved by setting a positive growth rate (inverse Landau damping) for one or several chosen wave numbers.To avoid the influence of boundary conditions and the finite size of the system, the growing mode wave number(s) should not be chosen either too small or too large.In the simulation runs analyzed in this article, we used a single pumping wave at approximately 1/5 of the Nyquist wave number.
In the initial stage of the evolution we observe a Langmuir wave spectrum consisting of several peaks, forming a cascade that transfers energy toward the low-k region.After a certain time, the modulational instability threshold is reached and a strong, low wave number component appears in the spectrum.At the same time solitary structures called cavitons (areas of density depletion and stronger electric field) are formed (see Fig. 1).Here the system reaches an "equilibrium" turbulent state, in which the total integral intensity |E| 2 dx remains approximately constant in time.It should be mentioned at this point that the dynamics of the cavitons in the 1-D case is fundamentally different from the physically, more realistic 3-D case.In three dimensions the cavitons are known to be unstable: as the modulation instability develops, they become deeper and narrower with time until they collapse (Zakharov, 1972).However, in the one-dimensional case the ponderomotive force is balanced by the pressure of the expelled plasma and stable cavitons are formed.
This way of simulating the Langmuir turbulence meets one additional difficulty, due to the finite size of the simulation box.When taking a very low damping rate for long-wave Langmuir and ion-sound waves (corresponding to physical conditions in space plasmas), one obtains significant growth of Langmuir and ion-sound waves in the low-k part of the spectrum.If the box size were very large, the longest wavemode in the system would be determined by the balance between characteristic time of the energy transfer toward the small wave numbers and the damping rate.In the case of a finite simulation box, the minimum wave number is defined by the size of the system, and this cutoff will influence the dynamics of the low-k part of the spectrum, where a significant amount of energy is cumulated, due to the weak dissipation of low wave number modes.
The cavitons in our simulations become very stationary and their spatial distribution in the simulation box is almost periodic (Fig. 2 shows their temporal evolution).This stationarity is in a good agreement with theoretical models (Rudakov, 1973) which describe the cavitons in the presence of ion-sound damping as forced non-propagating density fluctuations.This phenomenon was later confirmed by numeric simulations by Degtyarev et al. (1979), Al'terkop et al. (1976), Doolen et al. (1985) and many others.The quasiperiodic distribution of cavities is mostly a consequence of the finite size of the system and contributes to the excitation of standing waves, with wavelengths corresponding to the average spacing between the cavities.This effect contributes significantly to the overall dynamics of the system when the pumping growth rate is small.To ensure that this effect does not provide a major contribution to the evolution of the system, the pump-wave increment must be chosen to be large enough.
Our analysis is based on this stationary regime; we shall consider mainly one data set of 32 000 samples, sampled at the ion plasma frequency ω pi .The growth rate of the beam instability was set to γ = 0.01 ω pi at a single wave number of k = 0.066λ −1 D .The Landau damping rate was fixed to γ L k = Ck 5 and the constant C was chosen in such a way that the damping was effective at wave numbers k > 0.1 λ −1 D .The damping rate of ion-sound waves was set to a value of γ S k = 0.003 ω pi for all wave numbers.The simulation box size was 4000λ D and a 512-point FFT was used to perform the spatial Fourier transforms and to compute the convolutions.
Results obtained with other growth rates will be briefly discussed in Sect.7. Unless stated otherwise, we shall always subtract the time average of the variables prior to analysis.This average is essentially zero for the electric field, but is non-negligible for the density.
In Fig. 3 we depict the average power spectrum of the electric field and the density.We recognize the features described above, namely a peak at k = 0.066 in the electric field spectrum, which corresponds to the pump-wave (i.e. the energy input), and the next step of the cascade at k = −0.052.The presence of a small peak at k = 0.036 suggests another step, but this peak is hidden in the part of the spectrum where the modulational instability should be dominant.The two peaks in the density spectrum at k = ±0.118are signatures of ionsound waves that are produced by the Langmuir wave decay, as will be seen later.Our objective is to investigate the mutual interactions of these waves and quantify the nonlinear energy transfers among them.
Before we start with the wave analysis, it is appropriate to comment on the dispersion relation of the electric field waves shown in Fig. 4. The dispersion relation is inferred from the wave number-frequency power spectrum.It consists of two qualitatively different parts.In the positive frequency region the usual dispersion branch of Langmuir waves appears.In the negative frequency region most of the oscillation power is concentrated in a featureless region centered on the frequency ω = −0.02ω pi .These oscillations correspond mainly to non-propagating wave modes that are trapped inside the cavities.Note that since the electric field was demodulated (Eq.3), the positive (resp.negative) frequencies correspond to frequencies above (resp.below) ω pe .
This separation of the dispersion relation into a positive and a negative part is easily understood from the general form of the dispersion relation of the Langmuir waves In the presence of the slowly varying density fluctuations this relationship includes the frequency dependence upon local plasma frequency Since δn/n 0 < 0 inside the cavitons, the trapped wavemodes inside project to frequencies below ω pe and the propagating modes above ω pe .The negative frequency part is a manifestation of the trapped wave activity that cannot be described in the weak turbulence approximation and is closely related to the formation of the cavities.The frequency localization of these modes is determined by the depth of the cavities by virtue of Eq. ( 5).This result is consistent with a similar signature in ω − k space previously observed in Vlasov simulations by Goldman et al. (1996).This nonlinear frequency shift represents a bright example of strong turbulence characteristics that can be observed without special nonlinear analysis tools.

Higher-order spectral analysis
According to theory (i.e.Musher et al., 1995), for sufficiently strong damping of ion-sound waves and small amplitude of Langmuir waves, the low frequency density fluctuations are essentially forced by the electric field oscillations: Note that E is the field envelope, so the zero frequency in this plot corresponds to ω pe .The color scale is logarithmic.
The effect described by this equation is called the ponderomotive force.The properties of the Green functions G l,m determining the force are discussed in the cited works.If, in addition, we assume the validity of the weak turbulence approximation (i.e. if the characteristic time scales associated with the nonlinearity are much longer than those characteristic periods of the linear waves), then Eq. ( 1) can be approximated by Here, ω k is the complex frequency of the Langmuir waves, with its imaginary part corresponding to the growth rate or damping of the wave amplitude, and V kk 1 k 2 are quadratic coupling coefficients.This equation describes the evolution of high-frequency Langmuir oscillations in terms of nonlinear three-wave coupling of Langmuir and ion-sound waves.Specifically, it suggests that a dominant process in the current scenario is the decay of a Langmuir wave into a Langmuir and an ion-sound wave (l → l + s).By eliminating ρ k from Eq. ( 7) using Eq. ( 6) we can derive an equation describing the dynamics of the amplitudes of high frequency waves in a closed form: This equation makes it evident that there can appear phase relationships between different constituents of the wave spectrum.These relationships can be detected using an appropriate technique based on high-order spectra (HOS) which is (as shown by Kim and Powers, 1979) a relevant tool for studying such wave-wave interactions.
If, for example, an ensemble of waves with Fourier amplitudes U 1 , U 2 , . . ., U n interact along the resonance k 1 + k 2 + • • • + k n = 0 then, even though each single mode may have a randomly varying phase, there should exist a functional relationship between the phases of these waves.This property provides the basis for the use of HOS, to quantify the relative importance of possible interactions between Langmuir and ion-acoustic waves.
The mathematical theory of the HOS is based on the notion of cumulants, that are routinely used to describe the statistical properties of hydrodynamic turbulence (see, for instance, Monin and Yaglom, 1963) are defined as follows: Let u 1 (t 1 ), . . ., u n (t n ) be random functions and their characteristic function.Their joint cumulant is then defined by If the u j (t j ) are stationary in t j , the cumulants are functions of only n − 1 variables τ 1 , . . ., τ n−1 , where τ j = t j − t n .The cumulant spectrum F (k 1 , . . ., k n−1 ) is defined as the Fourier transform of C(τ 1 , . . ., τ n−1 ).As shown, in for example, Kim and Powers (1979), the cumulant spectra can be expressed as cumulants of the respective Fourier components of u j , matching the resonance condition In our analysis, we will make use of two important properties of cumulants which follow directly from the above definitions: 1.If an ensemble of random variables U 1 , . . ., U n can be divided into two statistically mutually independent groups, then the cumulant C[U 1 , . . ., U n ] is zero; 2. Cumulants (for n > 1) are invariant with respect to constant shifts in the variables For the purposes of our work we consider the two lowest order HOS, namely the bispectrum and the trispectrum, which are, respectively, defined as the second and third order cumulant spectra.Their application to wave-wave interactions is based on the property 1 mentioned above.If for example, the cumulant of three Fourier modes (the bispectrum), whose wave numbers (or frequencies) satisfy the resonance condition k 1 + k 2 = k 3 , is non-zero, the phases of the Fourier components are coupled to each other.In the same way, the trispectrum quantifies four-wave interactions.We must stress, however, that a non-zero HOS is necessary but not a sufficient condition for having nonlinear wave interactions (Pécseli and Trulsen, 1993).This point will be discussed later in this section.
For an ensemble of four complex Fourier modes (written as U k , V k , W k , X k ), the estimates of the bispectrum and trispectrum (Kim and Powers, 1979) can be written as where n = k + l − m and where we assume that The latter is possible without loss of generality due to property 2. The angular brackets denote averaging over any ensemble of statistical realizations of the random variable.In the following analysis, we will assume ergodicity of the process and replace the ensemble averaging by time averaging (Frisch, 1995).According to Eqs. ( 7) and ( 8), the relevant HOS for studying our system are the cross-bispectrum B EEρ and the autotrispectrum T EEEE .It is often more convenient to work with normalized quantities (respectively, bicoherence and tricoherence), in which the dependence on the spectral amplitude is eliminated, so that the absolute value allows one to quantify the relative power involved in the coupling (Kim and Powers, 1979).The bicoherence and tricoherence are respectively, defined as The absolute values are bounded between 0 (no correlation) and 1 (full correlation).Note that there exists other slightly different normalizations of the bi-and tricoherence (Kravtchenko-Berejnoi et al., 1995).
The modulus of the cross-bicoherence b EEρ is shown in Fig. 5.The clear peak at (k, l) = (0.066, −0.052) attests a strong phase coupling between the two strongest Langmuir modes l 0.066 , l −0.052 and s 0.118 .Since the wave k = 0.066 is the energy source in our system, it is likely that this phase coupling results from a Langmuir wave decay l 0.066 → l −0.052 + s 0.118 , wherein energy is transferred toward a lower wave number.The peak in the density spectrum at k = 0.118 corresponds to the ion-sound wave produced by this decay.In the same way, the peak at k = 0.036 is supposedly generated by the next step of the cascade.However the bicoherence does not allow one to resolve this decay, due to the relatively small amplitude of this wave compared to the surrounding wide band spectrum.
In Fig. 5 the diagonal region of relatively high bicoherence (reaching levels up to 0.3) corresponds to a phase coupling of the remaining part of the E k spectrum to the strong peaks in the low wave number part of density spectrum.This phase coherence is of a different nature.As mentioned before, in this region the threshold is reached for the modulational instability, which then becomes the prevalent process.From Eq. ( 6), which describes the effect of the ponderomotive force on the plasma density, it follows that this effect D , but only the part of the principal domain corresponding to significant power density is shown.Not all of the usual symmetries of bicoherence apply, due to the asymmetry of the power density |E k | 2 .The peak value of 0.45 at (0.066, -0.052) identifies a strong three-wave phase coupling of the type l 0.066 ↔ l −0.052 + s 0.118 .contributes to the bicoherence (Eq.15).This point will be elaborated upon in Sect. 5.
According to Eq. ( 8), the Langmuir wave decay can also be investigated using the auto-tricoherence t EEEE (k 1 , k 2 , k 3 ), which quantifies the phase coherence resulting from the four-wave processes of the type l + l → l + l. Figure 6 (left panel) shows the auto-tricoherence corresponding to the wave interaction involving the pump-wave l 0.066 + l (k 1 +k 2 −0.066) ↔ l k 1 + l k 2 .This figure reveals that the auto-tricoherence reaches significant values only if k 1 or k 2 is close to −0.052, which is the wave number of the second strongest Langmuir wave in the spectrum.The main difference with respect to the bicoherence is that we now observe a one-parametric set of phase coupled wave-modes l 0.066 + l k ↔ l −0.052 + l k+0.066−0.052, which is parametrized by the wave number k.
The right panel of Fig. 6 shows the auto-tricoherence for the case of k 3 = 0.064 (slight detuning with respect to the pump wave) corresponding to the process l k 1 +l k 2 ↔ l 0.064 + l (k 1 +k 2 −0.064) .No significant phase coupling is observed in this plot or in the remaining part of the three-dimensional tricoherence domain.This overall low level of tricoherence in regions that do not involve the pump wave again supports the validity of our description.
To summarize, the bi-and tricoherence analysis reveals the existence of a phase coupling between the strongest wave modes of the system.This result supports the weak turbulence approach to the description of the system, but at the same time several difficulties emerge.First, the phase relationship does not necessarily mean that there exists energy exchange between the modes, though it gives a strong argument in favor of that.Second, HOS do not allow one to The above threshold values indicate the presence of four-wave phase coupling corresponding to process l k 1 + l k 2 → l 0.066 + l (k 1 +k 2 −0.066) .In the right panel the tricoherence for fixed k 3 = 0.064λ D is depicted (a wave mode not included in the cascade).
The tricoherence plot shows no significant features and contains only random noise.
resolve any subsequent steps in the energy cascade (in our particular case, we cannot unambiguously say whether the wave at k = 0.036 is a product of Langmuir wave decay or not).These problems can be overcome by using Volterra modeling.

Volterra model analysis
Let us carry the HOS analysis one step further by introducing the concept of the energy transfers.We follow the computational framework developed by Ritz and Powers (1986) and later improved by Ritz et al. (1989).Since the electric field is expected to follow Eq. ( 7), we can describe its dynamics in terms of a general Volterra model and estimate the linear and quadratic kernels k and klm from the simulated data.These kernels contain all the pertinent information about the physical process.The field E k is expressed as (Ritz et al., 1989) and the time derivatives of |E k (t)| and φ k (t) which appear on the left side of Eq. ( 17) after substitution of Eq. ( 18) are approximated by finite differences.We then obtain an equation with Equation ( 19) can be viewed as a causal input-output model, which predicts the output electric field E(t + δt) from the inputs E(t) and ρ(t).The unknown Volterra kernels Q k lm and L k are estimated using a least-squares scheme (see the Appendix).We selected δt to be the sampling period (δt = 2π ω −1 pi ) and each pair of subsequent samples is used as one statistical realization.The number of unknown coefficients was restricted to cover only the part of the spectrum where the power density is significant (approximately −0.07 < k < 0.07 for the electric field and −0.13 < k < 0.13 for the electron density).Our experience shows that a larger range does not improve the estimates.
More physical insight into the nonlinear processes can be gained by introducing the kinetic equation for the spectral power P k = E * k E k , which is easily derived from Eq. ( 17) The energy transfer functions are the main quantities of interest, since they quantify the spectral power change that is due to nonlinear interactions.Energy transfer functions are more informative than HOS, since they reveal the magnitude of the energy transfer and, more importantly, its direction.Negative values of T k lm correspond to the decay l k → l m + s l , while positive values correspond to the inverse process l m + s l → l k .
In Fig. 7 we plot the energy transfer function T k1 k 2 ,k 1 −k 2 .This quantity can be directly interpreted as the rate of change of the Langmuir wave energy at k = k 1 , due to nonlinear interactions of the Langmuir wave at k = k 2 with the resonant ion-sound wave at k = k 1 − k 2 .Note the strong energy transfer from the pump-wave to the next step of the cascade l 0.066 → l −0.052 + s 0.118 (the peaks at (0.066, −0.052) and (−0.052, 0.066)).We now have direct evidence for an energy cascade toward smaller wave numbers.This result had been anticipated by HOS analysis, but only Volterra modeling attests it in an unambiguous way.Note also how the energy transfer function reveals the next step l −0.052 → l 0.036 + s −0.088 of the cascade (weaker peaks at (0.036, −0.052) and (−0.052, 0.036)).The bicoherence analysis was unable to properly resolve this.
The strong energy transfers that we observe in the central part of the plot cannot be given a reasonable physical interpretation.As explained in Sect. 2 the low wave number part of the spectrum is strongly influenced by the non-physical effects of the finite simulation box.Furthermore, the strong turbulence effects driving the behavior of the cavitons cannot be adequately explained in terms of wave-wave interactions and do not follow our model Eq. ( 17).Therefore, we cannot expect the Volterra model to properly resolve the wave-field properties.As we shall see shortly below, the large energy fluxes estimated by the model in that region are a consequence of the model inadequacy and the poorly posed inverse problem.

Interpretation of the results
In the following we summarize the physical significance of the results presented so far and use these results to draw some conclusions on the underlying dynamics of our system.From the power spectrum, the dispersion relation and the bicoherence results, we can conclude that two qualitatively different processes act on the dynamics of the system.One is the energy cascade, whose dynamics has been explained and quantitatively described by Volterra model analysis.This process is a well-known, weak turbulence effect.One of the necessary conditions that is required for Eq. ( 17) to hold is the validity of the random phase approximation (RPA) E k E * l ≈ δ k−l .To check the RPA, we plot in Fig. 8 the normalized autocorrelation coefficient In this plot we observe that if either k or l is close to 0.066 or −0.052, then the value of the correlation coefficient c(k, l) becomes close to zero for all k = l in an agreement with the RPA.The RPA can indeed be safely accepted for waves that are in the cascade range, thereby confirming the Volterra model approach for the decay instability.However, the validity of the RPA becomes questionable in the middle diagonal region of the figure, where |c(k, l)| reaches up to the level of 0.2.This is one of the reasons for the failure of the Volterra modeling in that region.
The deviation from the RPA essentially comes from the modulational instability, which is responsible for the formation of cavitons and the evolution of the low wave number part of the spectrum.The finite size of the simulation box and their consequences (described in Sect.2) also contribute to this effect.
Based on our analysis we can now compare the influence of the modulational and decay instabilities on the evolution of the system.The bicoherence (Fig. 5) shows the phase correlation patterns for both instabilities.The clear lines of low bicoherence (at k 2 = 0.066 and k 1 = −0.052) in Fig. 5 correspond to a weak coupling of the strong modes of the Langmuir cascade to the low wave number part of the spectrum.This signifies that the phase correlation resulting from ponderomotive forcing (as described by Eq. 6) is negligible for these modes and their dynamics is completely driven by the wave decay.The situation is different for the weaker wave at k = 0.036 where both effects contribute with comparable strength.Unlike the bicoherence, the Volterra model is able to reveal the weak decay energy transfer to this wave number.This separation of the wave-wave effect from the background strong turbulence phase coupling is achieved by explicitly assuming the form of interaction (Eq.17). 1  The low wave number phase coupling corresponding to Eq. ( 6) is, to a large extent, responsible for the failure of the model in this area.The technique for estimating the coupling coefficients is based entirely on the phase correlations (all the input parameters in the estimation procedure have a form of high-order spectra) and since these correlations are a result of effects not included in the model, we obtain invalid results.The ill-conditioning of our regression problem is partly responsible for the high magnitude of these erroneous energy transfers, as shall be seen in the following section.

Statistical validation
Higher order statistical quantities are known to be prone to errors.Therefore, there remains an important issue to validate the statistical significance of our results.We shall see that the results obtained for the dynamics of the Langmuir wave cascade are significant, but that in the low wave number region the Volterra model analysis is biased by the model's inadequacy.
1 Note that in the procedure of estimation of nonlinear coupling coefficients, fourth-order spectra are used (see the Appendix) which provide extra input information in addition to the bispectrum.Let us first consider finite sample size effects.HOS are very sensitive to the lack of statistics, especially if the magnitude of the coherence functions is much smaller than 1.The variance of the bicoherence (Eq.15) can be estimated (Kim and Powers, 1979) as where M is the number of independent statistical realizations.In our analysis we use 32 000 samples, but these are not independent.It is possible, however, to estimate the number of independent realizations from the ergodic theorem (Frisch, 1995), by introducing the integral time scale defined by T int k = ∞ 0 R(τ ) dτ , where R(τ ) is the autocorrelation function ( t means a sum over all time samples) For our data set T int k ranges from 40 to 100 sampling periods, depending on the wave number.We have, therefore, taken M = 32000/T int k = 320 in order to obtain an upper bound for the variance.In this case the standard deviation of the bicoherence reaches about 0.055.Thus, in areas where the bicoherence is significant (typically in excess of 0.3), the relative error is less than 15%.
The expression for the variance of the tricoherence (Kravtchenko-Berejnoi et al., 1995) has exactly the same form as for the bicoherence (Eq.24), where b is replaced by t.The relative error of the largest tricoherence peaks we observe is about 30%.We conclude that the peaks associated with the Langmuir wave cascade are statistically significant.
To gain more insight into the statistical significance of the power transfers, we tested the results against those obtained with surrogate data (Kantz and Schreiber, 1997).The purpose of this test is to determine if our results are really a consequence of a nonlinear deterministic process or if they could result from some linear phenomena (e.g.nonstationarity) that has not been properly taken into account.As has been said before, all the inputs of the regression procedure for estimating the energy transfers are high-order spectral moments.If we take into account the non-zero mean ρ k of the density fluctuation spectrum (giving ρ k = δρ k +ρ k ), then the Volterra model (Eq.17) can be rewritten as The last term on the right side corresponds to the nonlinear coupling, but the central one represents a non-resonant (k = m) linear process.Surrogate analysis was used to separate the two parts and compare their contribution to the energy transfer.We created 30 data sets derived from the original one by phase randomizing the density fluctuation (but otherwise keeping the power density structure and keeping the phase information of the electric field).By this we should have mostly eliminated the coupling corresponding to the last term of Eq. ( 26), while the middle one should stay unaffected.After that, we have computed the energy transfer functions for all of these data sets and carried out a simple statistical analysis of the resulting 30 realizations.It is not sufficient to use a single realization, because the resulting energy transfers depend on the particular randomization.Figures 9 and 10 show the mean energy transfers and the standard deviation obtained from this analysis.Comparing these plots to Fig. 7 we may check that in the range of the Langmuir cascade (the peaks in the top left and bottom right corners of Fig. 7), the mean energy transfers of the surrogate data reaches at most 15% of the original values.We conclude that the dynamics of the peaks is really associated with the nonlinear last term of Eq. ( 26).The central part of the plot, on the contrary, shows relatively large values, which confirms the negligible contribution of nonlinear wave interactions in that region.
Since the standard deviation of the energy transfer function can easily be estimated from our model, it is useful to consider its value as an additional test.In the Langmuir wave cascade range, the standard deviation is small, thereby confirming the validity of our Volterra modeling.In the low wave number region, however, the standard deviation is relatively large, which means that the energy transfers are strongly dependent on the particular realization.This result, together with the relatively high condition number2 (ranging from 900 to 2700 as a function of k) of the linear system that one has to solve to obtain the Volterra kernels, further confirms the inadequacy of the Volterra model in that region.
7 Dependence of the nonlinear energy transfers on the growth-rate of the instability Finally, we briefly show how the energy transfers depend on the growth-rate of the pump-wave instability.Figures 11 and  12 show the spectrum and the energy transfer functions of waves generated in a simulation with a growth-rate increased from 0.01ω pi to 0.03ω pi .The deviation from weak turbulence approximation is stronger than in the previous situation.As a consequence, the energy cascade should have fewer steps; this can indeed be observed in Fig. 11, where only one step appears.As is demonstrated in the power transfer plot (Fig. 12), the magnitude of the corresponding energy flow from the pump wave (peak in the top left corner) is significantly stronger than in the previous scenario.
On the contrary, if the growth rate is decreased from 0.01ω pi to 0.003ω pi , more steps of the cascade can be identified in the spectrum in Fig. 13, namely the peak at k = 0.036 is now much more pronounced.Figure 14 clearly shows the energy transfers between the three peaks, but their magnitude is weaker proportionally to the decrease in the energy input to the system.The two sharp peaks in the density spectrum correspond to stationary oscillations in the area between the cavities.Their presence is a simulation artifact that appears, due to the limited number of cavitons in the finite simulation box.The strong coupling of the cavitons to the Langmuir oscillations (in the sense of Eq. 6) introduces a significant error in the power transfers in Fig. 14.

Conclusions
The key result of this study is that HOS and Volterra modeling are appropriate statistical tools for gaining a better understanding of wave-wave processes in weak and in strong turbulence.These techniques allow one to identify the dominant wave-wave interactions and to analyze their properties.For the successful application of Volterra modeling, it is important that the model adequately describes the interactions of interest.However, this method allows one to extract the spectral energy transfers, even in conditions where other physical mechanisms affect the dynamics, as we demonstrated in the case of Langmuir turbulence, where both weak and strong turbulence effects were present.The analysis should, therefore, be always complemented by careful statistical validation, to distinguish the physically significant results from a bias introduced by the inadequacy of the model.
In this work we applied the techniques to simulation data, where we took advantage of unprecedented spatial resolution and of the large statistical content of the data set.The methods were also previously applied to experimental data (Kim and Powers, 1979;Ritz et al., 1989;Dudok de Wit et al., 1999), where the spatial resolution was limited to several observation points.The Volterra model analysis of experimental data requires the field to be measured at two or more spatial points (to estimate the spatial derivative), which is often a limiting factor, especially in the context of satellite observations of space plasmas.From this point of view, the CLUSTER experiment (involving four satellites) opens new perspectives for the application of similar models, possibly generalized to two or three dimensions.

Appendix -Linear dispersion and growth rate estimation
In this Appendix we describe the actual method used to estimate the linear and quadratic coupling coefficients and compare the growth rate and linear dispersion estimates with the true values.
The unknown Volterra kernels, Q k lm and L k , are estimated using the least-squares method by minimizing the error func- where Ê k is the output field as predicted by the model and E k is true experimental value of the output.Since the system (Eq.17) is linear in the unknowns, Q k lm and L k , the problem of minimizing (Eq.A1) is reduced to a multiple linear regression for N C unknown parameters Q k lm and L k .For each fixed k we obtain an over-determined set of N S − 1 (N S stands for the number of samples in the data set) equations where U k is a rectangular (N S − 1) × N C matrix.Each realization (one pair of subsequent samples) allows one to form one equation by substituting the experimental values of E k , E k , ρ k into Eq.( 19).This linear regression problem is solved using the conventional approach (Golub and Van Loan, 1989), where Eq. ( A2) is multiplied on the left by U * k , to obtain a linear system with a square matrix Note that the matrix U k contains the values of E k and ρ l E m for all the samples in our data set.Therefore, the coefficients of the matrix 1 N S −1 U * k U k are, in fact, the HOS of the types E k E * l , E k ρ l E * m and ρ n E k ρ l E * m .The estimation of the quadratic coupling coefficients is, therefore, entirely based on the HOS.The system (Eq.A3) with a positive definite Hermitian matrix can be cheaply solved by Choleski decomposition (Golub and Van Loan, 1989).Nevertheless, in the quadratic case the number of unknowns rarely exceeds several hundreds; therefore, its solution does not represent a significant numerical obstacle.The main numerical complication related to the solution of this problem is the illconditioning of the square matrix, as noted in Sect.6.
Once the coefficients L k and Q k lm are determined, it is possible to compute the more physically relevant quantities k and k lm using the relations (Eq.20).In these formulas the linear phase shift appears between the subsequent samples δφ k = ω k δt which needs to be estimated first.Ritz et al. (1989) suggest to estimate this quantity by a linear approximation (A4) The linear dispersion relation ω k estimated in this way is depicted in Fig. A1a (red line) in comparison with the true value (green line).It is evident that due to the strong contribution of the nonlinear terms, this estimate, based on a linear approximation, is inapplicable.This point is often overlooked: a nonlinear model is needed to fit the linear contribution.Fortunately, the model (Eq.19) allows one to determine the linear part of the dispersion relation directly as Im[ln L k ], which can be easily seen from the third relation in Eq. ( 20) by taking into account that δφ k = Im[ k ]δt.The blue line in Fig. A1a represents this approximation and confirms that, in the case under consideration, this method is particularly effective.The linear coefficient k is conventionally decomposed as Here ω k = δφ k /δt gives the estimate for the linear dispersion relation, as discussed above, and γ k is the average linear growth-rate of the wave with wave number k. Figure A1b shows the estimate for γ k obtained from our quadratic model, compared to the true value and the linear approximation γ lin k = (ln | E k E k / E k E k |)/δt.The true linear growth rate consists of a single peak at k = 0.066 representing the pump wave and is almost zero everywhere else, because the Landau damping is not effective in the low wave number region.Here the deviation of our estimate from the true values is more noticeable, especially in the higher wave number region.This has to do with the fact that since |γ k | |ω k |, it is more sensitive to statistical errors and contributions of processes not explained by the model.Nevertheless, Fig. A1 shows that this estimate is still much better than the linear one, because a large part of the nonlinear contributions to the growth rate is filtered out.

Fig. 1 .
Fig. 1.A time snapshot of the electric field intensity |E(x)| 2 and density fluctuation ρ(x) showing the cavities in the density and corresponding electric field peaks.

Fig. 2 .
Fig. 2. Excerpt of the spatio-temporal dynamics of the electron density.The cavities move very slowly and they conserve their almost periodic spatial distribution.

Fig. 3 .
Fig. 3. Wave number power spectrum of demodulated complex electric field (in dimensionless units) and of normalized density fluctuations.The Nyquist wave number is 0.4 λ −1 D .

Fig. 4 .
Fig. 4.Power spectral density P ω,k of the electric field vs. wave number and frequency, showing the Langmuir branch (ω > 0) and the non-propagating oscillations (ω < 0).Note that E is the field envelope, so the zero frequency in this plot corresponds to ω pe .The color scale is logarithmic.

Fig. 5 .
Fig. 5.The cross-bicoherence b EEρ (k 1 , k 2 ) between ρ end E.The Nyquist wave number is 0.4λ −1D , but only the part of the principal domain corresponding to significant power density is shown.Not all of the usual symmetries of bicoherence apply, due to the asymmetry of the power density |E k | 2 .The peak value of 0.45 at (0.066, -0.052) identifies a strong three-wave phase coupling of the type l 0.066 ↔ l −0.052 + s 0.118 .

Fig. 6 .
Fig.6.The auto-tricoherence t EEEE (k 1 , k 2 , k 3 ).Left panel shows the tricoherence computed for a fixed value of k 3 = 0.066λ D , corresponding to the strongest Langmuir mode in the spectrum.The above threshold values indicate the presence of four-wave phase coupling corresponding to process l k 1 + l k 2 → l 0.066 + l (k 1 +k 2 −0.066) .In the right panel the tricoherence for fixed k 3 = 0.064λ D is depicted (a wave mode not included in the cascade).The tricoherence plot shows no significant features and contains only random noise.

Fig. 7 .Fig. 8 .
Fig. 7.The energy transfer functions T k1 k 2 ,k 1 −k 2 quantifying the nonlinear energy transfer from a Langmuir wave at k 2 to Langmuir wave at k 1 .

Fig. 9 .
Fig. 9.The mean energy transfer functions obtained by surrogate analysis of 30 data sets with phase randomized density fluctuation.The axes are the same as that of Fig. 7.
Fig. 10.The standard deviation of the energy transfer functions obtained by surrogate analysis of 30 data sets with phase randomized density fluctuation.

Fig. 11 .
Fig. 11.The wavenumber power spectra of the electric field and density for increased pump-wave growth rate γ = 0.03ω pi .

Fig. 12 .
Fig. 12.The quadratic energy transfer functions for the increased growth rate γ = 0.03ω pi .

Fig. 13 .
Fig.13.The wavenumber power spectra of the electric field and density for decreased pump-wave growth rate γ = 0.003ω pi .

Fig. 14 .
Fig. 14.The quadratic energy transfer functions for the decreased growth rate γ = 0.003ω pi .
Fig. A1.Estimates of the linear growth-rate and linear dispersion relation as obtained by the linear and quadratically nonlinear model.