## Abstract

We consider the propagation of periodic waves with initially narrow spatial spectra in a Kerr medium. The set of constants of motion closely related to the order amplitudes is introduced. It is shown that the spatial spectrum remains uniformly narrow with propagation for both self-focusing and self-defocusing nonlinearity. In addition, for sufficiently weak nonlinearity, initially strong orders always remain strong. Thus the problem is shown to be essentially finite dimensional and well approximated by a proper set of coupled ordinary differential equations.

© 2004 Optical Society of America

## 1. Introduction

The nonlinear Schrödinger equation (NLS) describes the propagation of light in Kerr media, both in two spatial coordinates (photorefractive crystals and other nonlinear materials) and in one spatial and one temporal coordinate (fibers). The problem of periodic initial conditions has attracted attention recently because of possible applications in high-frequency fiber lines, nonlinear wave guiding, switching, and so on. For discussion we will use the language of two-dimensional (2D) spatial propagation, but the results can be easily reformulated for one temporal and one spatial coordinate as well.

There is a formal solution of the periodic NLS based on the analysis of the associated linear problem spectrum, which generally has an infinite number of bands. If the number of bands is finite, then complicated analytical solutions that are almost periodic in propagation coordinates can be built [1]. It is not clear how this procedure can be applied in practice to a problem with given initial conditions.

On the other hand, the propagation of periodic patterns in nonlinear Kerr-type materials often can be modeled by considering the interaction between a limited number of spectral (diffraction) orders [2]. Physically, the interaction is produced by the beam diffraction on the dynamic gratings that it writes in a material. For the important case when the period of the pattern is small, the Bragg–Kogelnik approximation is often sufficient for calculations. In this approximation, the formation of higher orders does not occur, order intensities are conserved, and the nonlinear interaction is reduced to modifications of phase velocities for interacting beams. Processes such as self-focusing or self-defocusing that lead to the enrichment of the spatial harmonic spectrum evidently cannot be explained in this approximation, but numerics demonstrates that taking into account a relatively small number of orders (e.g., 3–4 or more) often gives an adequate description for limited propagation length. In the paper by Thyagaraja [3] motivated by the earlier numerical results in the context of fluid dynamics it was demonstrated that for the kth-order intensity |*ψ*_{k}
|^{2}, |*ψ*_{k}
|^{2}<*C*
_{2}/|*k*|^{2} with propagation-length-independent constant *C*
_{2}. This estimation ensures only slow diminishing of spectrum tails and, as was shown with numerics in Ref. [4], strongly overestimates the number of orders necessary for adequate description.

By considering higher integrals of motion we refine the result of Thyagaraja. We demonstrate that for narrow initial spectrum |*ψ*_{k}
|^{2}<*C*_{n}
/|*k*|
^{n}
for any *n* with propagation-length-independent constants *C*_{n}
. For small nonlinearity (both positive and negative), we prove additionally that initially strong orders always remain strong upon propagation; i.e., only the small part of their initial intensity is transferred to higher orders.

The NLS for smooth periodic initial conditions is essentially finite dimensional. Only a limited number of harmonics can be taken into account to obtain a generally adequate description of the propagation. The trajectories predominantly lie on *N*-dimensional tori in a 2*N*-dimensional phase space, with a number *N* determined by initial conditions and nonlinearity strength. The rest of the harmonics produce a fine “jitter” in a phase space; they slightly modify they invariant tori form and characteristic angular frequencies. This remarkable feature is due to the special structure of the NLS, which allows for an infinite set of conservation laws.

The paper is organized as follows. In Section 2 we reproduce known results on the monodromy matrix and classic conservation laws necessary for what follows. In Section 3, the estimations for order intensities are made. In Section 4, a set of conserved quantities are given in terms of harmonic amplitudes. Section 5 presents numerical illustrations. Finally, Section 6 contains conclusions and discussion.

## 2. Monodromy matrix and conservation laws

Here we reproduce some known results on the monodromy matrix and conservation laws, which are necessary for the main part of the paper. This section closely follows Chap. 1 in Takhtajan and Faddeev [5], and we use those authors’ notation. The NLS is written in the following form:

Complex wave function *ψ* is x periodic, *ψ*(*x*+2*L*, *t*)=*ψ*(*x*, *t*). We attribute to *t* the role of propagation spatial coordinate. Here *ϰ* is positive for the self-defocusing situation and negative for the self-focusing.

The associated linear problem is written as

with a two-component vector function *F*. Here *λ* is a spectral parameter, and the matrix *U* is

Here *$\overline{\psi}$* is complex conjugate. For the square root, the positive value is taken for *ϰ*>0 and the value with the positive imaginary part for *ϰ*<0.

The monodromy matrix *T*_{L}
is introduced as the matrix that gives the associated problem solution change over the *x* period so that *F*(*L*,*t*,*λ*)=*T*_{L}
(*t*,*λ*)*F*(-*L*,*t*,*λ*).

The monodromy matrix *T*_{L}
can be formally found with a solution of the Fredholm-type integral equation

where

Then *T*_{L}
(*λ*)=*T*(*L*,-*L*,*λ*), and we omit *t* dependence.

It can be shown that *tr*(*T*_{L}
(*λ*)) for any *λ* is a conserved quantity; i.e., if the *t* evolution of *ψ* is described by the NLS, the monodromy matrix trace remains the same. This permits us to find a set of conserved quantities by calculating the trace for the large real spectral parameter, namely,

Because *Tr*(*T*_{L}
) is a conserved quantity, the real coefficients *I*_{n}
are conserved upon propagation. These quantities were first found in Ref. [6] and constitute the classic set of integrals of motion for the NLS. They can be expressed as integrals of *ψ*,*$\overline{\psi}$* and their *x* derivatives, namely,

with

We need some of these integrals to be written explicitly. Some of their terms are modified from the expressions above with integration by parts (see also Ref. [6]). The intensity *I*
_{1} is

and we introduce the average intensity *I*.

The Hamiltonian *I*
_{3} is

The fifth integral is

## 3. Estimations for order intensities

Our aim is to establish limits for the transversal spatial coordinate spectrum of the wave function *ψ*

where *k* is integer and *K*=*π*/*L* is a wave vector.

For positive *ϰ* corresponding to self-defocusing, a simple estimation is possible. Both terms in *I*
_{3} [Eq. (12)] are then positive, which means that

thus

The same argument works for *I*
_{5}, giving the estimation with *k*
^{-4}. Nevertheless, it fails for *I*
_{7} where the term proportional to *ϰ* can be negative for some functions.

There is a more elaborate procedure that permits *t*-independent estimation that does not rely on the nonlinearity sign and gives a decrease faster than any power of *k*. We demonstrate this procedure for *I*
_{3} first. Let us introduce

The point *x*
_{0} is chosen to correspond to the minimum of |*ψ*(*x*)|^{2}. Then |*ψ*(*x*
_{0})|^{2}≤*I* and 0≤6 Φ(*x*)≤2*LI* for *x*
_{0}≤*x*≤*x*
_{0}+2*L*. The second integral in Eq. (12) can be then written as

The first term in Eq. (18) does not exceed 2*LI*
^{2}. To estimate the absolute value of the next two terms, we use the Cauchy–Schwarz inequality:

Then

Here we have used the upper limit for Φ and the intensity conservation. The same estimation is obtained for the third term in Eq. (18). Now it is seen that the upper limit for the first term in *I*
_{3} exists, because if the first term is large, the second is proportional to the square root of it. The uniform estimation can be obtained as follows:

Thus, for strong self-focusing nonlinearity, taking into account only the leading term in the Eq. (21),

where the dimensionless nonlinearity parameter *ν*=ϰ*I*(2*L*)^{2} describes the nonlinearity strength. This means that the typical spectrum width is proportional here to |*ν*|. For strong self-defocusing nonlinearity, the width, as follows from Eq. (16), is proportional to √|*ν*|. Thus the self-focusing nonlinearity produces more-pronounced widening of the spectrum in comparison with the self-defocusing one. This is to be expected, since for the self-focusing a modulation instability exists, which is absent otherwise. Equation (21) was first obtained in Ref. [3].

The left-hand side of the Eq. (18), as follows from Eqs. (20) and (21), is limited by *Q*
_{1}=2*LI*
^{2}+2(2*LI*)^{3/2}
*A*
_{1}. Then, similar to the derivation for *I*
_{3}, from the *I*
_{5} conservation it follows that

By induction one can show that all terms for higher integrals are limited as well. The proof for *I*
_{5} and the procedure for higher integrals are given in Appendix A.

## 4. Order-related integrals of motion

The results of Section 3 demonstrate that there is a uniform estimation for the order amplitudes with a decrease for large *k* faster than any negative power of |*k*|. This indicates that there is an exponential estimation |*ψ*_{k}
|^{2}<*C*exp(-*α*|*k*|) with constants *C* and *α* determined by the initial conditions and *ϰ*. By estimating the width of the exponential distribution, and comparing it with Eqs. (16) and (22) for large *ν*, it is seen that *α* scales approximately as |*ν*|^{-1} for self-focusing nonlinearity, and as *ν*
^{-1/2} for self-defocusing.

Thus, for a given nonlinearity strength and narrow initial spectrum, only a finite number *N* of lower Fourier harmonics is important, whatever the propagation length. The higher orders introduce uniformly small perturbation. Because the NLS is known to be completely integrable, the solutions are approximately restricted to *N*-dim tori in 2*N*-dim phase space of *N* complex amplitudes. The form of the tori is determined by conservation laws. One can use *N* first integrals of Zakharov–Shabat to this purpose. Nevertheless, we will introduce another set of conservation laws, which have simpler structure for small nonlinearity. They are obtained with standard iterations of the Fredholm-type integral equation [Eq. (4)] for an infinite discrete set of special values of *λ*, *λ*_{k}
=*kπ*/*L*.

When we introduce

the monodromy matrix trace is

with

and so on.

With the function *f*(*z*)=$\underset{-L}{\overset{z}{\int}}$
*dz′µ*_{k}
(*z′*), and integration of Eq. (26) by parts, is easy to show that

Thus for small nonlinearity the conserved quantities are simply order intensities, which corresponds to linear propagation. The *G*
_{1}+*G*
_{2}+… can be considered as a modulus square of generalized *k*-th Fourier harmonic, which is conserved for the nonlinear case. Different from the linear case, *G*
_{1}+*G*
_{2} +… do not diminish exponentially with *k*→∞ for smooth initial data. By comparison with the asymptotical expression (7), it is seen that this sum behaves as -*ϰ*
^{2}
${I}_{1}^{2}$
*L*
^{2}/2(*kπ*)^{2}. To obtain the set that is exponentially diminishing for *k*→∞, one can subtract the terms of Eq. (7) that behave as 1/*k*^{n}
for all n. After this, the term *O*(|*k*|^{-∞}) still remains, because the series in Eq. (7) are asymptotic, and not exact. The example of constant function *ψ*(*x*)=*const* for which the monodromy matrix can be found analytically demonstrates that this procedure can work only for sufficiently large *k*; otherwise, the series in Eq. (7) do not converge.

If expressed in Fourier harmonics, the *G*
_{2} includes combinations of the type *$\overline{\psi}$ _{k}*

_{1}

*ψ*

_{k}

_{2}

*$\overline{\psi}$*

_{k}_{3}

*ψ*

_{k}

_{4}, which are resonant; i.e., to obtain nonzero contribution, it is necessary that -

*a*

_{1}

*k*

_{1}+

*a*

_{2}

*k*

_{2}-

*a*

_{3}

*k*

_{3}+

*a*

_{4}

*k*

_{4}=

*a*

_{5}

*k*with

*a*

_{1}, …,

*a*

_{4}=0,1 and integer

*a*

_{5}. The explicit equation for this coefficient is complicated because there are many possible resonant combinations of different types.

The convergence of series Eq. (25) follows from the theory of integral equations. The usual proof [5] includes estimation with the maximal value of |*ψ*|. It is possible to estimate the terms with the Cauchy–Schwarz inequality as well (see Appendix B). In this case, one has

This proves the convergence, as does the usual procedure, but the estimation is made with the average intensity, which, different from the maximal value of |*ψ*|, is propagation independent.

From the estimation for *G*
_{2} and Eq. (28) it follows that if initially

the *k*th order preserves its identity, because the sum of the subsequent terms for *J*_{k}
will be smaller than *G*
_{1}. Thus, if initially the order is strong, and the nonlinearity, as determined by the *ν* parameter, is small, the amplitude evolves in a characteristic ring in a phase plane (see Fig. 2 for a numerical illustration). If the nonlinearity parameter *ν* is large, the initially strong orders lose their identity, and the ring becomes a circle. In this case there is still a finite number of strong orders, but this number is greater than the number of orders that were strong initially.

## 5. Numerical illustrations

To determine the actual form of the spectrum when the number of harmonics in it is initially limited, we have performed numerical calculations of propagation for initial conditions given by 2–4 lower Fourier harmonics with randomly taken phases and amplitudes.

The set of ordinary differential equations was obtained from the NLS by representing *ψ* as a Fourier series [Eq. (14)]. After substituting the series into the NLS and forcing all higher harmonics to be zero, the set of *N* differential equations in *N* variables is obtained, as, e.g., in Ref. [2]. For more efficient numerical solution we compensated for the linear part of propagation by the ansatz:

After this, the equations for *ψ′*_{k}
were solved numerically with the fourth order Runge–Kutta method. The consistency was checked by adding new orders to calculation and doubling the number of steps. The algorithm permits us to calculate with sufficient exactitude the solution with 20 orders taken into account for lengths up to 10^{4} or more. To investigate the diminishing of intensity with the number of diffraction orders, we plotted the maximal value along the propagation path for all harmonics involved. Different propagation lengths were considered. For the initial conditions investigated, the plots look quite similar. The typical result is presented in Fig. 1. The plots exhibit exponentially diminishing tails with a slope depending on the initial conditions and nonlinearity strength. As a function of propagation length there is a pronounced saturation. In fact, the points for propagation lengths *l*=10^{3} and 10^{4} are indistinguishable on the graph. The difference in logarithms for the highest harmonic is 0.004. The chosen value for self-focusing nonlinearity is already in a region of modulation instability.

For the nonlinearity parameter (4<|*ν*|<40), which is not very small, the slope of exponential diminishing in Fig. 1 scales well as |*ν*|^{-1}. For large n the number of orders that must be accounted for grows, which makes calculation more difficult. For negative nonlinearity of the same strength, the slope is larger than for positive one, and it diminishes more or less proportionally to |*ν*|^{-1/2}, although the correspondence is worse than in the self-focusing case.

In Fig. 2 we show the complex plane evolution for the strong order at a moderate nonlinearity. It is seen that the order intensity upon propagation varies in certain limits only, in accordance with the theory presented above. Equation (30) seems to underestimate the value of *ν* necessary for order destruction ~1 order of magnitude. It is possible, though, that other configurations of beams lead to closer values.

## 6. Discussion and conclusions

The NLS with periodic initial conditions is “essentially finite dimensional.” This means that the intensity concentrated in higher diffraction orders remains uniformly small upon propagation. If initial intensity is concentrated in a small number of orders *M*, the enrichment of spectrum related to self-focusing/defocusing occurs, but beyond this there is no further spreading. We show that the spectrum tails rapidly (exponentially) diminish. Thus the propagation can be modeled for any length by taking into account only the finite number *N* of diffraction orders. If nonlinearity is small, *N* is not much greater than *M*; for strong nonlinearity, especially the self-focusing type, the two can differ substantially. The harmonics beyond *N* introduce a uniformly small perturbation. Thus an adequate description can be obtained with a finite number of ordinary differential equations. If one uses an integrable finite-dimensional approximation (for example, by rewriting the integrable Ablowitz–Ladik set [7] for discrete Fourier harmonics), the solutions will lie on invariant *N*-dimensional tori in 2*N* dim phase space. The KAMtheorem implies that the distortions resulting from the rest of the harmonics will slightly deform these tori and introduce changes in phase velocities.

For rapidly diminishing initial conditions, the solitons and diffracting part are formed upon propagation. Our analysis in this paper shows that for the periodic initial conditions the bound structure of Fourier harmonics is formed around strong orders. The number of orders involved in this structure depends on the nonlinearity strength, which is essentially finite. This structure is characterized by *N* propagation constants giving characteristic phase velocities. If all or most of these frequencies are equal, the soliton-like structures similar to cnoidal waves are formed. For moderate nonlinearity the strong waves do not lose their identity; thus one can speak of the “dressed plane wave.”

Because only a finite number of characteristic frequencies exist, the total behavior is nearly recurrent. The recurrence was observed in a number of numerical experiments. This is also consistent with the fact that quasi-periodic finite band solutions form a dense set in a phase space. Our treatment gives a complementary picture, showing that the solution with smooth initial data is close to a quasi-periodic one. It is important that the number of independent frequencies can be estimated from the initial data, either with Eq. (21) or with higher integrals.

## Appendix A

The estimation of terms entering *I*
_{5} can be made as follows. Let us introduce Φ_{1}(*x*)=$\underset{x0}{\overset{x}{\int}}$|*ψ*|4*dz* with *x*
_{0} corresponding to the minimum of |*ψ*(*x*)|^{2}. The function Φ_{1}(*x*) is limited with *t*-independent constant *Q*
_{1}for *x*
_{0}≤*x*≤*x*
_{0}+2*L*, as we have already established considering the *I*
_{3} integral. Then an equation similar to Eq. (18) can be written for the final term in the Eq. (13):

The argument then proceeds exactly as in Eqs. (19)–(20), with the limits *Q*
_{1} and Φ_{1} instead of 2*LI* and Φ. Note that now we have in the equation corresponding to Eq. (20) the integral that we have already estimated for *I*
_{3}.

If we introduce Φ_{2}(*x*)=$\underset{x0}{\overset{x}{\int}}$
|*ψ*_{z}
|2*dz* with the same *x*
_{0}, the argument of the previous paragraph can be repeated for the two terms in Eq. (13) that have the first *x* derivatives. Again, we have a *t*-independent upper limit for Φ_{2} established at the earlier step for *I*
_{3} and can proceed along the lines of Eqs. (19)–(20).

Because all terms, except for the first in Eq. (13), are limited, and their sums are conserved, there is a t-independent limit for the first term as well. This gives *C*
_{4}/*k*
^{4} bound for |*ψ*_{k}
|^{2}.

In general, the proof by induction can be made as follows: The terms in *P*
_{2n-1}, Eq. (10), have the structure

with *φ* meaning either *ψ* or *$\overline{\psi}$*, and the number in parentheses is the derivative order. It is seen that *p* is even and that the sum *p*+*i*1+…+*ip*=2*n* for all terms belonging to *P*
_{2n-1}. We will call 2*n* the rank of a term. One can always obtain the higher derivative in terms of rank 2*n* smaller or equal than *n*-1 integrating the term by parts. We will call the term estimated if we have established that the integral of its modulus is limited upon propagation. Let us assume that all the terms of rank 2*n*-2 or smaller are already estimated, and consider the terms with rank 2*n*. The only term having two *φ*, after integration by parts, is *ψ*
^{(n-1)}
*$\overline{\psi}$*
^{(n-1)}. This is the one that we can estimate only after all the rest of the same rank are estimated. Consider the terms with *p*=4,6, … The construction of Φ is made in the following way. First, let us separate the term in pairs with even rank for each pair (the rank of the pair is sum of its derivative orders plus two). This can always be done, since the number of *φ* with even derivative is even (otherwise, the rank of the term will be odd). Then choose the pair with the smallest rank. The rank of this pair *r* is 2≤*r*≤*n*. The rest has the rank smaller or equal to 2*n*-2, and it is estimated by hypothesis. So one can build with it a Φ function similar to that of Eq. (17). The pair has a rank smaller than 2*n*, and, consequently, it is also estimated. Then one proceeds as in Eqs. (18)–(20). The derivative increases the pair’s rank by 1, so in the worst case the rank of 1 *φ* in the integrals of the equation equivalent to Eq. (20) is *n*, and the rank of another is 1. In this case the estimation gives the square root of the term ∫*ψ*
^{(n-1)}
*$\overline{\psi}$*
^{(n-1)} . All other possibilities permit estimation with terms of a smaller rank. Thus the estimation of all terms can be done by induction. This ends the proof.

## Appendix B

Equation (29) is obtained as follows. Let us estimate *G*
_{2} using the Cauchy–Schwarz inequality, Eq. (19). The square of the modulus of the final integral is

Then

and so on, finally until we obtain Eq. (29).

## Acknowledgments

The financial support of CONACyT project U39681 is gratefully acknowledged. We thank one of the Reviewers for providing the reference to early research on the subject [3, 4].

## References and links

**1. **Y. C. Ma and M. J. Ablowitz, “The periodic cubic Schrödinger equation,” Stud. Appl. Math. **65**, 113–158 (1981).

**2. **N. Korneev, A. Apolinar-Iribe, and J. J. Sanchéz -Mondragon, “Theory of multiple beam interaction in photore-fractive media,” J. Opt. Soc. Am. B **16**, 580–586 (1999). [CrossRef]

**3. **A. Thyagaraja “Recurrent motions in certain continuum dynamical systems,” Phys. Fluids **22**, 2093–2096 (1979). [CrossRef]

**4. **D. U. Martin and H. C. Yuen “Spreading of energy in solutions of the nonlinear Schrödinger equation,” Phys. Fluids **23**, 1269–1271 (1980). [CrossRef]

**5. **L. A. Takhtajan and L. D. Faddeev*Hamiltonian Methods in the Theory of Solitons* (Springer-Verlag, 1987).

**6. **V. E. Zakharov and A. B. Shabat, “Exact theory of two-dimensional self-focusing and one-dimensional self-modulation of waves in nonlinear media,” Sov. Phys. JETP **61**, 62–69 (1972).

**7. **M. J. Ablowitz and J. F. Ladik, “Nonlinear differential-difference equations and Fourier analysis,” J. Math. Phys. **17**, 1011–1018 (1976). [CrossRef]