# Breakdown of Modulational Approximations in Nonlinear Wave Interaction

## I Introduction

Modulational instability of high-frequency nonlinear waves is a common process in a variety of circumstances involving wave propagation in continuous systems. Modulational processes can be seen to occur in a wide range of physical situations, from nonlinear waves in plasmas  to nonlinear electromagnetic waves propagating in optical fibers . What usually happens in all those cases is that due to generic nonlinear interactions, the amplitude of a high-frequency carrier develops slow modulations in space and time. If the modulations are indeed much slower than the high-frequencies involved, one can obtain simplified equations describing the dynamics of the slowly varying amplitudes solely, the amplitude equations . In the present analysis we consider systems that become integrable in this modulational limit, a feature often displayed. Should this be indeed the case, no spatiotemporal chaos would be observed there. The basic interest then would be to see what happens when the approximations leading to modulational approximations cease to be satisfied. The paper is organized as follows: in §2 we introduce our model equation and discuss how and when it can be approximated by appropriate amplitude equations; in §3 we investigate the modulational process from the point of view of nonlinear dynamics; in §4 we perform full spatiotemporal simulations and compare the results with those obtained in §3; and in §5 we conclude the work.

## II Model equation, modulational approximations, and amplitude equations

In the present paper we focus attention on a nonlinear variant of the Klein-Gordon equation (NKGE) to investigate the breakdown of modulational approximations in the context of nonlinear wave fields. The NKGE used here reads

$$_t^2A(x,t)_x^2A(x,t)+\frac{\mathrm{\Phi }}{A}=0,$$
(1)   
($`_t/t,_x/x`$) where we write the generalized nonlinear potential as

$$\mathrm{\Phi }(A)=\omega ^2\frac{A(x,t)^2}{2}\frac{A(x,t)^4}{4}+\frac{A(x,t)^6}{6},$$
(2)   
$`\omega `$ playing the role of a linear frequency which will set the fast time scale. The remaining coefficients on the right-hand-side were chosen to allow for modulational instability and saturation; we shall see that while the negative sign of the second fulfills the condition for modulational instability, the positive sign of the third provides saturation. The choice of their numerical values is arbitrary, but our results are nevertheless generic. The NKGE is known to describe wave propagation in nonlinear media and the idea here is to see how the dynamics changes as a function of the parameters of the theory: wave amplitude, and time and length scales.

Let us first derive the conditions for slow modulations. We start by supposing that the field $`A(x,t)`$ be expressed in the form

$$A(x,t)=\stackrel{~}{A}(x,t)e^{i\omega t}+\mathrm{complex}\mathrm{conjugate}.$$
(3)   
Then, if one assumes slow modulations and discards terms like $`_t^2\stackrel{~}{A}`$ and the highest-order power of the potential $`\mathrm{\Phi }`$, one obtains

$$2i\omega _t\stackrel{~}{A}(x,t)_x^2\stackrel{~}{A}(x,t)3|\stackrel{~}{A}(x,t)|^2\stackrel{~}{A}(x,t)=0,$$
(4)   
which is, apart from some rescalings, the Nonlinear Schrödinger Equation - we shall refer to it as NLSE here - an integrable equation.

This is no novelty; it is known that the modulational approximation is obtainable when there is a great disparity between the time scales of the high-frequency $`\omega `$ and the modulational frequency, we call it $`\mathrm{\Omega }`$, such that terms of order $`\mathrm{\Omega }^2\stackrel{~}{A}`$, when compared to $`\omega ^2\stackrel{~}{A}`$, can be dropped from the governing equation. The magnitude of the modulational frequency can be estimated as follows. Consider Eq. (4) and suppose a democratic balance among the magnitude of its various terms; $`\omega _t\stackrel{~}{A}_x^2\stackrel{~}{A}\stackrel{~}{A}^3`$. Then one obtains for $`_t\mathrm{\Omega }`$,

$$\frac{\mathrm{\Omega }}{\omega }(\frac{\stackrel{~}{A}}{\omega })^2.$$
(5)   
It is thus clear that the modulational approximation is valid only when $`\stackrel{~}{A}\omega `$, since this condition slows down the modulational process causing $`\mathrm{\Omega }\omega `$. The next question would be on what is to be expected when the modulational approach ceases to be valid. Before proceeding along this line, let us mention that a stability analysis can be performed on Eq. (4). One perturbs an homogeneous self-sustained state with small fluctuations of a given wavevector $`k>0`$ (we choose $`k>0`$ here, but the theory is invariant when $`kk`$) and after some algebra one concludes that: (i) the field in the homogeneous state, let us call this field $`A_h`$, is given by

$$A_h=a_oe^{i\frac{3a_o^2}{2\omega }t}$$
(6)   
where $`a_o`$ is an arbitrary amplitude parameter and where the exponential term should be seen as providing a small nonlinear correction to the linear frequency $`\omega `$; and (ii) the perturbation is unstable when

$$k<k_{tr}\sqrt{6}a_o.$$
(7)   
with maximum growth rate at

$$k_{max}\frac{k_{tr}}{\sqrt{2}}.$$
(8)   
When unstable, the homogeneous state typically evolves towards a state populated by regular structures which can be formed precisely because the underlying governing NLSE, Eq. (4), is of the integrable type as mentioned earlier.

## III Beyond the modulational approximation

To advance the analysis beyond the modulational regimes we start from the basic equation, Eq. (1), but do not use approximation $`_t^2\stackrel{~}{A}\omega ^2\stackrel{~}{A}`$ leading to Eq. (4) and eventually to condition $`\stackrel{~}{A}\omega `$. The idea is precisely to examine what happens as the ratio $`\stackrel{~}{A}/\omega `$ of Eq. (5) grows from values much smaller, up to values comparable to the unit.

Our first task is to examine how the regular and well known modulational instability analyzed in the previous section comes directly from Eq. (1). To do that, let us write a truncated solution as the sum of an homogeneous term plus fluctuations with wavevector $`k`$, $`A=A_h(t)+A_1(t)(e^{ikx}+e^{ikx})`$, $`A_h`$ and $`A_1`$ real. The truncation, that discards higher harmonics, is legitimate within linear regimes, but fortunately we shall see that it is not as restrictive as it might appear even in nonlinear regimes. The general idea favouring truncation here is that in the reasonable situation where modes with the fastest growth rates are more strongly excited, relation (8) indicates that second harmonics are already outside the instability band since $`2k_{max}=\sqrt{2}k_{tr}>k_{tr}`$. Under these circumstance one would be led to think that the most important modes would be the homogeneous and those at the fundamental spatial harmonic. We shall actually see that the truncation provides a nice and representative approach to the case of regular regimes.

Note that we are already considering the amplitudes of the exponential functions as equal. This results from a simplified symmetrical choice of initial conditions and is totally consistent with the real character of Eq. (1). After some lengthy algebra, one finds out that the coupled nonlinear dynamics of the fields $`A_h`$ and $`A_1`$ is governed by the Hamiltonian

$$H=\frac{p_h^2}{2}+\frac{\omega ^2q_h^2}{2}\frac{q_h^4}{2}+\frac{2q_h^6}{3}+\frac{p_1^2}{2}+\frac{\chi ^2q_1^2}{2}\frac{3q_1^4}{4}+\frac{5q_1^6}{3}$$

$$3q_h^2q_1^2+10q_h^4q_1^2+15q_h^2q_1^4,$$
(9)   
where $`\chi ^2\omega ^2+k^2`$, $`q_h=A_h/\sqrt{2}`$, $`q_1=A_1`$, and where the $`p`$’s denote the two momenta conjugate to the respective $`q`$-coordinates.

The Hamiltonian (9) can be informative. As a first instance it can be used to determine the stability properties of the homogeneous pump, as mentioned before. To see this, assume that in average $`q_hq_1`$ and solve the dynamics perturbatively. In zeroth order, one would have the following Hamiltonian governing the $`(p_h,q_h)`$ dynamics:

$$h_o=\frac{p_h^2}{2}+\frac{\omega ^2q_h^2}{2}\frac{q_h^4}{2}+\frac{2q_h^6}{3}.$$
(10)   
With help of action-angle variables ($`\rho ,\mathrm{\Theta }`$) for the zeroth-order part and conventional perturbative techniques the solution reads

$$q_h=\sqrt{\frac{2\rho }{\omega }}\mathrm{cos}(\omega t\frac{3}{2}\frac{\rho }{\omega ^2}t),$$
(11)   
if $`\rho `$, the amplitude parameter, is not too large. Note that the oscillatory frequency undergoes a small nonlinear correction which is determined by the quartic term in $`q_h`$ of the Hamiltonian $`h_o`$.

Next we consider the driven Hamiltonian controlling the dynamics of the canonical pair $`(p_1,q_1)`$

$$h_1=\frac{p_1^2}{2}+\frac{\chi ^2q_1^2}{2}+3q_h^2q_1^2,$$
(12)   
where we recall that the pair $`(p_1,q_1)`$ describes the inhomogeneity of the system. One again introduces action-angle variables $`(I,\theta )`$ to rewrite the linear Hamiltonian (12) in the form

$$h_1=\chi I+12\frac{\rho }{\omega }\frac{I}{\chi }\mathrm{cos}^2\theta \mathrm{cos}^2(\omega t\frac{3}{2}\frac{\rho }{\omega ^2}t),$$
(13)   
from which we obtain the resonant form

$$h_{1,r}=\frac{1}{2\omega ^2}(k^23\rho )I+\frac{3\rho I}{2\omega ^2}\mathrm{cos}(2\varphi ),$$
(14)   
with $`\varphi =\theta (\omega 3\rho /2\omega ^2)t`$, and where use is made of the approximation $`\chi \omega +k^2/(2\omega )`$ valid when $`k\omega `$. Now consider the stability of a small perturbation $`I0`$. For this initial condition $`h_{1,r}0`$. Since $`h_{1,r}`$ is a constant of motion, solutions for arbitrarily large values of $`I`$, what would indicate instability, are possible only when $`|\mathrm{cos}(2\varphi )|1`$. From the resonant Hamiltonian (14), this demands

$$k<\sqrt{\frac{6\rho }{\omega }},$$
(15)   
and also indicates that maximum growth rate for $`I`$ occurs at $`k_{max}=\sqrt{3\rho /\omega }`$. Comparisons of temporal dependence of Eqs. (3), (6) and (11) shows that $`\rho /\omega =a_o^2`$, and that conditions (7) and (15) are therefore one and the same as they should be. In other words, starting from the full nonlinear wave equation, we recover the typical results naturally yielded by the NLSE.

But the Hamiltonian (9) gives further information because apart the truncation, it is not an adiabatic approximation like Eq. (4); it can therefore tell us whether the reduced dynamics, if unstable, is of the regular or chaotic type. The interest on this issue comes from the fact that the reduced dynamics usually helps to determine the spatiotemporal patterns of the full system: while regular reduced dynamics is associated with regular structures - frequently a collection of solitons or soliton-like structures, chaotic reduced dynamics is associated with spatiotemporal chaos. The correlation has its roots on the so called stochastic pump model . The model states that intense chaos in a low-dimensional subsystem of a multidimensional environment can make the subsystem act like a thermal source, irreversibly delivering energy into others degrees-of-freedom. In the limit of a predominantly regular dynamics undergoing in the reduced system, irreversibility is greatly reduced and energy tends to remain confined within the subsystem. On the other hand, in the limit of deeply chaotic dynamics with no periodicity at all, energy flow out of the subsystem is fast and there is not even much sense in defining the subsystem as an approximately isolated entity. The intermediary cases are those more amenable to a description in terms of the stochastic pump . We point out, however, that in all cases, even in the chaotic one, analysis of the reduced subsystem serves as an orientation on what to be expected in the full spatiotemporal dynamics. In our case, the appropriate subsystem is precisely the one we have been using. This is so because it is the smallest subsystem comprising the most important ingredients of the dynamics: the homogeneous and the only linearly unstable modes.

We now examine these points in some detail. Let us consider $`\omega =0.1`$. When $`A_{h,o}\omega `$, $`A_{h,o}A_h(t=0)`$, the modulational approximation can be used to determine the instability range, $`k<k_{tr}`$. A surface of section plot based on the two-degrees-of-freedom Hamiltonian (9) produces Fig. (1), where we record the values of $`(q_1,p_1)`$ whenever $`p_h=0`$ with $`dp_h/dt>0`$, and where we take $`q_1=0.01q_hq_h`$ and $`p_h=p_1=0`$ to determine the total unique energy of the various initial conditions we launch in the simulations. Note that the particular “seed” initial condition introduced above is included in the ensemble of initial conditions launched and represents small perturbations to an homogeneous background.

In Fig. (1a) we take $`A_{h,o}/\omega =0.5`$ and $`k=2k_{max}>k_{tr}`$. One lies outside the instability range and the figure reveals that the origin, which represents the purely homogeneous state $`q_1=A_1=0`$, is indeed stable. The remaining panels are all made for $`k=k_{max}`$ and increasing values of the ratio $`A_{h,o}/\omega `$. One sees that for such values of $`k`$ and $`A_{h,o}`$ not only the origin is rendered unstable, as it also becomes progressively surrounded by chaotic activity. The inner chaotic trajectory issuing from the origin - i.e., the trajectory representing the modulational instability - is encircled by invariant curves, but the last panel already shows that some external chaos, here represented by scattered points around the invariant curves, is also present. At some amplitude $`A_{h,o}=A_{critical}`$ slightly larger than the one used in panel (d), invariant curves are completely destroyed. Above the critical field, inner orbits are no longer restricted to move within confined regions of phase-space - these orbits are in fact engulfed by the external chaos seen in panel (d). External chaos here is a result of the hyperbolic point positioned at the local maximum of the generalized potential $`\mathrm{\Phi }(A)`$, at $`A_10.1`$. One gross determination of the critical field based on simple numerical observation yields $`A_{critical}/\omega 1/1.62`$, although the corresponding transition from order to disorder in the full simulations may not be so sharply defined. Yet, one may expect that in both low-dimensional and full simulations, the transition should occur at comparable amplitudes.

## IV Full spatiotemporal simulations

To make the appropriate comparisons, we now look into the full simulation of Eq. (1). Full simulations are made through a discretization of the spatial domain via a finite difference method. The dynamics is resolved temporally by means of a sympletic integrator, as was the purely temporal Hamiltonian. The results are quite robust and energy is conserved up to $`10^6`$ parts in one. In Fig. (2) we display the spacetime history of the quantity $`Q=\sqrt{\omega ^2A(x,t)^2+\dot{A}(x,t)^2}`$ \- we plot this quantity because it becomes a constant of motion in the limit where we discard nonlinearities and inhomogeneities. Initial conditions are the same as the seed initial condition used in Fig. (1), but the control parameters differ. In the case of Fig. (2a) where we choose $`A_{h,o}/\omega =0.1`$ so as to safely satisfy $`A_{h,o}A_{critical}`$ and $`A_{h,o}\omega `$, it is seen that the spatiotemporal dynamics is very regular. The homogeneous state is unstable, but only periodic spatiotemporal spikes can be devised. This is the regular spatiotemporal dynamics so typical of the integrable NLSE and this kind of dynamics agrees very well with the low-dimensional predictions of the reduced Hamiltonian (9). The fact that only one single structure can be found along the spatial axis at any given time, indicates that the dynamics is singly periodic (period-1) along this axis and thus basically understood in terms of the reduced number of active modes (homogeneous plus fundamental harmonics) of Hamiltonian (9). We now move into the vicinity of the critical amplitude $`A_{critical}`$ discussed earlier. Under such conditions, one may expect to see the effects of spatiotemporal chaos. The value $`A_{h,o}/\omega =1/1.625`$ is chosen in Fig. (2b), where full simulations indeed display a highly disorganized state after a short regular transient. Regularly interspersed spikes can no longer be seen and spatial and temporal periodicities are lost, which characterizes spatiotemporal chaos. In this regime many modes become active (we shall return to this point later) and the reduced Hamiltonian (9) fails to provide an accurate description of the full dynamics. However it still provides a good estimate on the point of transition. A more throughful examination of the transition in the full simulations suggests that the critical field there is a bit smaller - a value close to $`\omega /1.95`$; for smaller values we have not observed noticeable signals of spatiotemporal chaos even for much longer runs than those presented here. The much larger oscillations executed by $`Q`$ in chaotic cases (see the legend of Fig. (2)) is a direct result of the destruction of the invariant curves seen in Fig. (1). In the absence of invariant curves, initial conditions are no longer restricted to move near the origin.

It is thus seen that there are limits to an integrable modulational description of the dynamics of a wave field. The limits are essentially set by the parameter $`A_{h,o}/\omega `$. If it is much smaller than the unit, the modulational description is valid and one can expect to see a collection of spatiotemporal periodic structures being formed as asymptotic states of the dynamics. On the other hand, as the parameter approaches the unit, nonintegrable features are likely to be seen. In particular, regularity does not survive for very long, and fluctuations with various length scales appear in the system. This is a regime of spatiotemporal chaos which fundamentally involves the presence of nonlinear resonances between the frequency of the carrier, $`\omega `$, and the intrinsic nonlinear modulational frequency, $`\mathrm{\Omega }`$.

Due to the presence of chaos, one is suggested that the transition involves an irreversible energy flow out of the reduced subsystem. If one computes the average number of active modes

$$<N^2>\frac{\underset{n}{}n^2|A_n|^2}{_n|A_n|^2},$$
(16)   
where the amplitudes are defined in the form

$$A_n=\underset{j}{}A(x_j,t)e^{inkx_j},$$
(17)   
one obtains the plots shown in Fig. (3).

In Eq. (17), “$`j`$” is the discretization index. In Fig. (3a) we use the same conditions as in Fig. (2a). This panel shows that in rough terms, energy keeps periodically migrating between the homogeneous mode (when $`\sqrt{<N^2>}0`$) and the fundamental harmonic (when $`\sqrt{<N^2>}1`$). Conditions of Fig. (3b) are the same as those of Fig. (2b); one sees that as the ratio $`A_{h,o}/\omega `$ grows, periodicity is lost, and that energy flow out of the initial subsystem into other modes becomes clearly irreversible. In Fig. (3c) we use the same previous conditions with exception of $`A_{h,o}`$ which we take $`A_{h,o}=\omega /1.95`$. This slightly smaller, but not too small, value of the initial amplitude allows to observe the slow diffusive transit of energy during the initial stages of the corresponding simulations. During this stage one can actually look at the subsystem as an energy source adiabatically delivering energy into other modes - the concept of the stochastic pump applies more appropriately in those situations.

## V Final conclusions

To summarize, in this paper we have studied the breakdown of modulational approximations in nonlinear wave interactions. We have analyzed a nonlinear Klein-Gordon equation to draw the following conclusions. Adiabatic or modulational approximations are accurate while the high-frequency of the carrier wave keeps much larger than the modulational frequency. Under these circumstances the full spatiotemporal patterns are regular as is the dynamics in the reduced subsystem where the system energy is initially injected. There is no net flow of energy out of the reduced subsystem into the remaining modes.

On the other hand, when both frequencies become of the same order of magnitude, the reduced subsystem undergoes a transition to chaos. Correspondingly, the spatiotemporal patterns of the full system become highly disordered and energy spreads out over many modes. The correlation between the low-dimensional and high-dimensional spatiotemporal chaos has its roots on the stochastic pump model . According to the model, a low-dimensional chaotic subsystem may act like a thermal source, delivering energy in a irreversible fashion to other degrees-of-freedom of the entire system. Spectral simulations performed here indicates that this seems to be the case with the present setting.

The transition to chaos involves a noticeable increase in terms of wave amplitude. This takes place when the invariant curves of Fig. (1) are destroyed, allowing for the merge of the external and internal chaotic bands into one extended chaotic sea. This merging of chaotic bands is actually a result of reconnections involving the manifolds of the unstable fixed point at the origin, and other hyperbolic points associated with the curvature of the generalized potential $`\mathrm{\Phi }`$ . When both manifolds reconnect, inner trajectories issuing from the origin start to execute the large and irregular oscillations seen in Fig. (2). More detailed studied of the process is under current analysis.

## ACKNOWLEDGMENTS

This work was partially supported by Financiadora de Estudos e Projetos (FINEP), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), and Fundação da Universidade Federal do Paraná - FUNPAR, Brazil. S.R. Lopes wishes to express his thanks for the hospitality at the Plasma Research Laboratory, University of Maryland. Part of the numerical work was performed on the CRAY Y-MP2E at the Supercomputing Center of the Universidade Federal do Rio Grande do Sul.