Branched flow and caustics in nonlinear waves

Rogue waves, i.e. high amplitude fluctuations in random wave fields, have been studied in several contexts, ranging from optics via acoustics to the propagation of ocean waves. Scattering by disorder, like current fields and wind fluctuations in the ocean, as well as nonlinearities in the wave equations provide widely studied mechanisms for their creation. However, the interaction of these mechanisms is largely unexplored. Hence, we study wave propagation under the concurrent influence of geometrical (disorder) and nonlinear focusing in the (current-modified) nonlinear Schrödinger equation. We show how nonlinearity shifts the onset distance of geometrical (disorder) focusing and alters the peak intensities of the fluctuations. We find an intricate interplay of both mechanisms that is reflected in the observation of optimal ratios of nonlinearity and disorder strength for the generation of rogue waves.


Introduction
Rogue waves are extreme fluctuations in random wave fields with much higher amplitude than the average wave height. Their occurrence is an expression of the fact that the amplitude or intensity statistics of the waves in the field follow a heavy tailed distribution. Probably the most sensational manifestation are giant ocean freak waves of more than 20 m height, which have been confirmed to exist by measurements for the first time only little more than 20 years ago [1]. Heavy tailed intensity distributions and thus rogue waves, however, can be found in a variety of different contexts including optics [1] and microwave transmissions [2,3], as well as the propagation of sound [4] and tsunami waves [5].
Several physical mechanisms have been identified that can lead to those heavy-tailed intensity distributions [6], two of which are arguably the most prominent ones: on the one hand this is the modulation instability occurring in nonlinear wave equations with focusing nonlinearity. In this context, rogue waves are supposed to be manifestations of special, spatially and temporally localized breather solutions of the wave equation [7,8]. On the other hand it is the formation of branched flows which has been found in the propagation of linear waves in weakly scattering but correlated random media and leads to strong spatial focusing. Branched flows are closely connected to the occurrence of random caustics [9,10], i.e.singularities in the density of the rays approximating the wave field(see e.g. [11]).
Breathers have been extensively studied both theoretically and experimentally in optics and water waves [12][13][14][15][16][17]. Likewise have branched flows been studied in a variety of systems ranging from transport of electrons in the two-dimensional electron gas scattered by the disorder potential of impurities [18][19][20], via the scattering of wind driven ocean waves by ocean currents [21,22] to the focusing of tsunamis by small variations of the ocean floor topography [5]. Not only do branched flows in general lead to heavy-tailed intensity distributions [10,23], but experiments on microwave transmissions through disordered arrays of scatterers have shown that in combination with fluctuating sources they can lead to breather-like isolated rogue events [2].
Despite this progress a full quantitative theory of the statistics of rogue waves, especially on the oceans, is still missing. One important ingredient needed to achieve this, we believe, will be a thorough understanding of the interplay of linear and nonlinear focusing mechanisms, which are typically both present in the ocean wave dynamics. Especially the formation of branched flows in the presence of modulation instability needs to be understood. First numerical observations have indicated that the combination of the two effects can indeed lead to higher rogue wave probabilities than the individual mechanisms [22,24].
With this article, we want to contribute to the systematic study of branched flows in nonlinear wave propagation, using the nonlinear Schrödinger equation (NLS) in the presence of weak disorder as the paradigmatic model equation. Figure 1 shows an illustrative example of a branched flow created by an initially plane wave propagating over a weak disorder potential, comparing the linear to the nonlinear wave dynamics in the same potential landscape. The precise model we use is described in section 2. Since in two-dimensional wave flows the most fundamental structure is the cusp caustic, we will then give a very brief review of branched flows in linear waves and its close connection to the formation of random caustics in section 3. Because of this close connection we start our analysis of the nonlinear wave propagation by studying the impact of nonlinearity on a single cusp caustic created by a curved wave front in the absence of any potential in section 4, complementing our numerical findings by perturbative nonlinear ray theory. In sections 5 and 6 we will then study the propagation of nonlinear waves in weak disorder potentials.

Model
In this article we will study the (2 + 1) dimensional NLS [25] A t where A is the (complex) wave envelope and the potential term Ṽ is a weakly fluctuating correlated random field representing e.g.a turbulent current field in the ocean or a correlated random, weakly refractive medium. This equation models the modulation of weakly nonlinear waves in different systems [1], including the propagation of water waves in the deep ocean, where it is denoted as the current-modified NLS. The random field V k V 0 = has zero mean and variance V 2 á ñ which we express as V k 2 The parameter ò, which we will use in the remainder of this article, measures the relative disorder strength, compared to the kinetic energy of an incoming (linear) plane wave with wave vector k 0 . The random potential is further characterized by a spatial correlation function which we choose to be Gaussian with width ℓ c , i.e.
For example in the ocean, the parameter  is given by the ratio of the root mean square velocity of the background current and the speed of the incoming wave, which is on the order of 10 −3 [22]. The dimensionless nonlinearity g˜is given by a K 2 0 0 2 ( ) , where a K 0 0 is the wave steepness of the carrier wave with typical amplitude a 0 and wavenumber K 0 . For deep water ocean waves the steepness can grow up to a few times 10 −1 [1].
We assume that the amplitude of the random potential is small, which allows us to use the paraxial approximation Figure 1. Intensity for two realizations of a branched flow for g=0 (bottom panel) and g=−0.1 (top panel) with k 0 =20π in the same potential V with ò=0.007, ℓ c =1. The distance to the first caustic, i.e. to the maximum of the averaged scintillation index å, is shown by the golden dashed line, as well as 10 times this distance. The nonlinearity has a strong influence on the branched flow, i.e. the onset distance is decreased through the self-focusing nonlinearity and the structures are sharpened.
In our simulations we solve equation (3) using a Fourier split-operator method following [26] with a typical resolution Δx=Δyλ 0 /10 (λ 0 =2π/k 0 ) and periodic boundary conditions in y-direction. Our main goal is to study the interplay of random and nonlinear focusing. To this aim, we study initially plane waves propagating in the random potential landscape. For a few examples, we qualitatively compared the solutions obtained in the paraxial approximation with much more costly simulations of the non-equilibrium steady states observable in the full 2 1 d + ( ) equation (1) in the presence of a planar wave source and absorbing boundary conditions using the methods described in [26]. In all cases we found very good agreement. An example is shown in figure 2. In contrast to optics, a plane wave initial condition is not very realistic in the ocean. Here, one rather expects a superposition of plane waves where the propagation direction and the wavenumber vary in a certain range. This will affect the probability distribution of the wave intensities, and in general reduces the chance to observe rogue waves, as described in [21,22,24] for geometrical and nonlinear focusing mechanisms individually. However, in order to focus our study on the principles of the interplay of caustics and branched flows with nonlinear focusing mechanisms, we limit the initial conditions to a single plane wave.
As mentioned above, branched flows are tightly linked to the formation of random caustics [9,10], i.e. singularities in the ray description of the wave field. We, therefore, will first study how individual caustics are affected by a (small) nonlinearity in the wave equation before we turn to branched flows. Caustics are characterized in catastrophe optics by their fundamental normal forms [11], which parameterize the variation of the optical path length in the vicinity of the caustic. In two-dimensional wave fields, only two types of caustics can occur: fold and cusp caustics. In a cusp point two fold caustics are spawned, therefore we can restrict our discussion to the optical cusp catastrophe. To create a single cusp caustic we need an initial condition u 0 (y)=u(0, y) in equation (3), which in the paraxial approximation captures the curvature of the wave front that leads to the desired optical path length variation given by the normal form of the cusp. An appropriate initial condition, which corresponds to a plane wave that has just passed through a suitable phase screen, is given by 3  The parameters r, a and β jointly control the shape of the cusp and the distance of the cusp point from the phase screen as will become apparent in section 4.1. The maximal phase difference along the initial condition is βk 0 . Here, we keep α=1000 and β=2π fixed.
The free, linear propagation (i.e. V(x, y)≡0 and g=0) of the initial condition(4) leads to focusing into the desired cusp structure in the intensity u 2 | | as can be seen in the center panel of the left column of figure 3 (this effect we will denote as geometrical focusing in the following). The upper and lower panels show the impact of weak self-focusing (g<0) and self-defocusing nonlinearities (g>0) on the cusp. The right column of figure 3 illustrates that equivalent structures can be seen in branched flows.
As might be expected, we can observe that in both cases a self-focusing nonlinearity leads to higher peak intensities, sharper structures, and earlier focusing whereas a self-defocusing nonlinearity leads to decreased peak intensities with blurred out structures and delayed focusing compared to the linear case. This effects can be seen even more clearly in figure 4 which shows cuts along the central horizontal line of the cusp structures of figure 3. The maximum intensity clearly increases for g<0 and decreases for g>0 and the distance to the first maximum decreases and increases respectively.
Before we study the impact of the nonlinearity on geometrical focusing and branched flows quantitatively in the later sections we will first very briefly review some statistical properties of branched flows in linear waves in the next section.

Scaling of branched flows (linear waves)
Branched flows show universal behavior for a wide range of random media, i.e.fundamental properties of the flow are independent of the details of the dynamics and the random medium and only depend on very few characteristic parameters like the variance 2  and the (integral) correlation lengths ℓ c of the fluctuations in the medium [5,23,28]. Similar to diffusion in random media which is characterized by the mean free path, branched flows scale with a characteristic length scale that is proportional to the mean distance to the first caustic x 0 c . In isotropic media 4 this length scales as which for small ò is much shorter than the mean free path that scales like ℓ c /ò 2 . How can one measure the branching length in a wave simulation or an experiment? Caustics are singularities in the ray field which approximates the wave front propagation, and their position is well defined and easily measured in numerics. In wave fields, however, they appear as more or less diffuse regions of increased intensity. In the case of a single cusp caustic created by a curved wavefront, we can use the position of the absolute maximum of the intensity to define the relevant length scale. In the complex wave fields of branched flows, as the ones of figure 1, we need to take another way. We have shown earlier [3,5] that the average fluctuation strength of the wave field (averaged over many realizations) quantified by its relative variance, the scintillation index (see section 5), is a useful measure. As a function of propagation distance from the source, in a branched flow the scintillation index shows a peak. The distance of this peak from the source defines a suitable length scale that we will use as x 0 c in the following. In this peak region the strongest branches are visible in the individual realizations.

Geometrical focusing
Caustics are well studied and described in geometric and linear wave optics for rays and linear waves, whereas there are only few studies about geometrical focusing in the NLS or water waves, see e.g. [27,30]. In this section, we study how a single cusp caustic, generated by a curved wavefront as described in section 2, shifts with the nonlinearity in the wave equation as illustrated in figures 3 and 4. More specifically, we want to quantify and explain how the distance to the maximum varies. This variation is shown in figure 5 for different curvatures of the initial wave front as a function of the nonlinearity g. For self-focusing nonlinearities the distance decreases. For defocusing nonlinearities it first increases but then surprisingly decreases again. To gain a deeper understanding of this unexpected decrease for positive g and the functional form for the expected decrease for negative g we investigate the dynamic in a semiclassical approximation. . Intensity for a cut through y=0 for different g with r=0.8 for the wave fields shown in figure 3. The onset distance is decreased and the intensity is increased for self-focusing nonlinearities, whereas the onset distance is increased and the intensity decreased for defocusing ones.

Semiclassical approximations
In the linear dynamics the rays propagate from the curved wave front through free space and reach the cusp at a distance x c To assess the effect of the nonlinearity in a perturbative manner, we study the propagation of rays (or trajectories) that are following Newtons equations of motion in a potential that is given by the solution A 0 of the linear wave equation times the nonlinearity, i.e.
where the initial conditions y y 0 and y 0 i ( ) have to be chosen such that the rays are perpendicular to the initial wave front, i.e.
in the paraxial approximation. Figure 6 shows examples of the resulting ray fields for different g. A clear qualitative difference can be seen between focusing and defocusing nonlinearity. For focusing nonlinearity the behavior does not qualitatively change with g: a caustic is found along the central line. While for a small defocusing nonlinearity the rays still form a caustic on the central line, for stronger defocusing nonlinearity the caustic splits into two, which lie off the central line. Along the central line the rays are initially focused, but before they can cross and form the caustic, the curvature of the wave front changes and the rays are defocused. Only later two smaller cusps form left and right of the central line. The numerical solutions of the nonlinear wave propagation show, however, that the wave intensity maximum is still on the central line. We therefore suppose that the position where the rays become parallel (and thus the wave front starts to defocus) is a good approximation for the distance to the maximum of the wave field. In the ray dynamics the cusp point or the parallel point respectively, depending on g, can be found numerically using the stability matrix. We compare them to the position of the maximum in the actual nonlinear wave field in figure 7. The curves show rough qualitative agreement and follow the same trend, but do not match quantitatively. However, we will show that the turning points, i.e.the positive nonlinearities at which the distance to the maximum starts to decrease again can be well predicted. Furthermore, we will show that for self-focusing nonlinearities the shape of the curve can well be approximated analytically.
A solution for the cusp is known and is given by the Pearcey integral [11]. Using this expression the perturbed distance to the caustic can in principle be calculated analyzing the stability matrix of the central trajectory. The latter is simply a straight line (traversed with varying speed), but using the full form of the integral, however, the resulting equations still appear not to be solvable analytically. Therefore, we approximate the curvature of A 0 2 | | (which enters the expression of the stability matrix, see appendix A) by a constant 0  < , which simplifies the equations sufficiently. The constant ñ implicitly depends on the initial wavefront, i.e.the parameters r, α and β, and will be used as a fit parameter. This allows us to determine a simple functional dependence for x c g as a function of the strength of the nonlinearity g. As shown in appendix A, the resulting estimate for the distance to the caustic is Figure 5. Distance to the maximum of A 2 | | for k 0 =2π together with a fit (lines) to the perturbatively found functional form equation (8). The vertical lines indicate the distance at which rays start to separate in the perturbative ray dynamics, given by g b = | | . Numerical results are presented by the symbols.
Using the same approach the distance to the parallel point can also be calculated and for g>0 we find In this approach we made two major approximations. First the distance was calculated based on ray dynamics and secondly the expression A y 2 0 2 ¶ | | along the central ray was approximated by a constant. Note that already in the linear case the wave intensity maximum is not at the caustic position but beyond, therefore a constant offset was included in the fit. Thus here and later, when we analyze branched flows, we use a fitting function of the form

Branched flow
We will now turn our attention away form the single caustic created by a curved wave front to the branched flow field created by the diffraction of waves by correlated disorder. In the latter an ever increasing number of random caustics is created.

Scintillation index
As mentioned above, a helpful quantity to characterize the wave flow is the scintillation index. It measures the relative strength of the intensity fluctuations and is defined as ) | | is the wave intensity. Here the averages are taken over realizations of the disorder. If the initial condition is a plane wave propagating in x direction, the scintillation index Σ=Σ(x) is a function of x only, i.e. a function of the propagation distance. Typical curves are shown in figure 8, for three different disorder strengths. In the numerics the ensemble average is taken from 1024 realizations of the random potential (and over the y position).
In the linear case (g=0) the scintillation index initially increases strongly and reaches a maximum when the strongest branches occur, then it slowly decays to one at large distances 5 . A scintillation index of Σ=1 is what you expect for a random superposition of plane waves, i.e. when the intensity is exponentially distributed [31,32].
Let us now examine how the scintillation index curves change with nonlinearity. In this and the subsequent section we fix k 0 =20π, and ℓ c =1. We see in figure 8 that the maximum of the scintillation index, i.e.the magnitude of the strongest intensity fluctuations, varies with the nonlinearity. It increases for self-focusing and decreases for defocusing ones, consistently for the different disorder strengths (a)-(c). But there is also a qualitative difference in the shape of the curves. After the initial peak due to the onset of branching, a second wider peak develops with increasing self-focusing nonlinearity. The extend of this second peak, however, strongly depends on the disorder strength. In (a) for the strongest disorder the scintillation index is still approaching one at large distances, signaling an exponential intensity distribution as in the linear case. But in the two topmost curves one can clearly see a second maximum. This behavior is more pronounced for smaller disorder strength in (b), where the decay to 1 is much slower. In the top curve the second maximum becomes an extended plateau, with nearly the same magnitude as the first peak. In (c) the second increase of the scintillation index is strongest and the magnitude rises above that of the initial peak. Here, even at the end of the simulation region, the intensity distribution will decay much slower than exponential. In (c) the linear case g=0 does not show a fully developed branched flow, as can be seen by the absence of a peak in the scintillation index. The nonlinearity, however, intensifies the fluctuations and sharp caustics appear. This is in accordance with the caustic sharpening observed in [33], which studied caustics created in a light-field passing through a phase screen and then propagated through a (homogeneous) nonlinear medium. For very weak disorder the self-focusing nonlinearity increases the intensity fluctuations and recovers the peak that was washed out in the linear case. In contrast, for the intermediate disorder strength in figure 8(b) we see that defocusing nonlinearities apparently leads to the disappearance of caustics. This is similar to the behavior for a geometrically focused wave, in which the defocusing prevents the crossing of the rays, at least on the central line.
It is apparent that the influence of the nonlinearity on the shape of the scintillation index is stronger for smaller disorder strength ò. Furthermore, one can see that the distance to the maximum depends on the nonlinearity. We will quantify this dependence in section 5.3. The qualitative differences can also be seen in the examples of branched flow fields of figure 1. The focusing nonlinearity leads to sharper branches of higher intensity. They are much more stable and persists over longer distances. Also the branches can interact and reflect one another. These kind of interactions are known for nonlinear waves and are a prominent feature of soliton-solutions. This suggests that the structures present in the branched flow are some kind of realization of these solutions. Different studies from optics have compared peaks appearing due to randomly perturbed phases of the initial wave and found that the shapes match known solutions quite well [12,13], though their complexity is higher [34].
In the next two subsections we want to quantify the changes in the scintillation index induced by nonlinearity.

Fluctuation strength
Let us first study the change in fluctuation strength as characterized by the maximum Σ max of the scintillation index. The lower panel of figure 9 shows Σ max as a function of g for the same three disorder strength as in figure 8. That the maximum in the scintillation index is indeed the consequence of the formations of caustics is strikingly confirmed in figure 9 which shows the dependence of the maximal intensity A max 2 | | of caustics generated by curved wave fronts as a function of g for three appropriately chosen initial curvatures of the wave front parameterized by r (see equation (4)) in the upper panel and the maximum of the scintillation index for different disorder strength in the lower panel. The correspondence of the two cases is very clear. The maxima decrease for positive g and initially increase for g<0 but then a turning point can be observed. For small curvatures and correspondingly small ò one can clearly see that the maxima decrease again. The same behavior can also be seen for larger ò only for stronger negative nonlinearities g<0 outside the shown simulation range. This means that depending on the disorder strength there is an optimal g in the sense, that the fluctuations are most intense. This is an interesting finding and shows that the interplay of the geometrical and nonlinear focusing is highly nontrivial. We will strengthen this point of view by studying intensity and rogue wave statistics in section 6. But first we want to conclude our examination of the scintillation index by studying the onset distance of the branching behavior.

Onset distance
As in the linear case, we define the onset distance of the branched flow as the distance x c from the source to the maximum of the scintillation index. The change of this distance with nonlinearity is shown in figure 10. Here we can see, qualitatively, the same behavior as for the single cusp caustic. For self-focusing nonlinearities the distance decreases and it initially increases for defocusing ones. In contrast to the geometrical focusing, we do not observe a decrease of the distance for stronger defocusing nonlinearities, but instead the peak in the scintillation index disappears, as can be seen in figure 8(b). This originates from the disappearance of the caustics. This is underlined by figure 3, in which for g>0 one can see, that the high intensity regions are smoothed out.
Again we argue that the scintillation maximum is caused by the formation of caustics and therefore we expect that the distance x c follows the same functional g-dependence as that of the individual cusp. In figure 10 we therefore fit expression equation (8) for self-focusing nonlinearities to the simulation results. In addition we mark for which g>0, from the results of section 4.1 and the parameters of the fit, we expect the single cusp to disappear. Overall we find a good correspondence.
So far we lack a proper understanding of the second peak in the scintillation index. One would assume a connection to modulation instability, however, its position varies with disorder strength, suggesting an interplay of the two effects.

Intensity and rogue wave statistics
Let us finally come back to the most intriguing question that was already raised in the introduction: how does the interplay of random and nonlinear focusing influence the rogue wave probability?
When linear waves propagate in weakly diffracting, correlated random media the formation of branched flows leads to heavy-tailed intensity distributions and therefore to a strong increase of the probability of rogue waves compared to the expectations from a random wave model with exponential intensity distribution [2,10,21,23]. A scintillation index larger than one is a consequence of this heavy-tailed distributions. At large distances from the source we have seen in section 5 that the scintillation index is approaching one from above and indeed the intensity probability density function (pdf) is approaching an exponential [32]. Therefore in a branched flow of linear waves there is typically only a certain range of propagation distances (near the peak distance x 0 c that we have discussed earlier) when rogue wave probabilities are strongly increased. This behavior of the linear waves is confirmed in the right column of figure 11 for the strongest disorder strength (and can also be found for intermediate disorder strength; not shown). The figure shows the intensity pdf compared to an exponential decay at two different distances from the source: at the scintillation peak x c and 10 times this distance. While the pdf for g=0 is clearly heavy tailed at x c it becomes exponential at 10x c . But more over the figure shows the influence of the nonlinearity: positive g lead to a faster decay of the pdf, while negative g lead to more heavy-tailed distributions and heavy-tailed intensity distributions prevail for much longer propagation distances.
The left column of figure 11 shows the same curves but for the smallest disorder we studied, where the linear branched flow is not fully developed (as described in section 5.2) and does not show a heavy-tailed but an exponential intensity pdf at both distances. However, focusing nonlinearities restore the heavy-tailed distribution.
To make the influence of the nonlinearity easier to grasp quantitatively, we calculated the probabilities of freak waves (i.e. rogue waves). We used a simplified definition and waves of intensity larger than I I 4.4 2 ) [21], in which I 0 á ñ is the mean intensity for g=0, are named freak waves. The probability P(I>I f ) as a function of the nonlinearity is shown in figure 12. In addition we present P I I* > ( ), with I I 4.4 0 * = á ñin figure 12. We show, as above, curves for the three different disorder strengths and at two different propagation distances, i.e.at the scintillation peak and 10 times that far. In accordance to what we saw for the maximal intensity of the curved wavefront and the maximum of the scintillation index there appears to be an optimal nonlinearity for the different disorder strengths, where the rogue wave probabilities are highest, clearly observed for the smallest ò for P(I>I f ) at distance 10x c and for P I I* > ( ). If there is a maximum for the larger ò but at stronger nonlinearities, outside of our simulation regime, we can only suspect. It is very clear, however, that the interplay of random and nonlinear focusing in the formation of rogue waves is highly non-trivial.

Conclusion and outlook
In conclusion, we have studied the impact of self-defocusing and self-focusing nonlinearities in the wave dynamics on single cusp caustics as well as branched flows in diffracting random media in the realm of the NLS using numerical and perturbative methods. We analyzed the impact on the onset position of high intensity fluctuations and their intensity statistics.
The most interesting findings concern the intensity statistics and the formation of rogue waves in branched flows in the presence of self-focusing nonlinearity. For small nonlinearities it is clear that the region of strong branching, characterized by a pronounced peak in the scintillation index, is the main realm for rogue wave formation. It would thus appear that caustic formation enhanced by nonlinear focusing is the prime origin of rogue waves. But the nonlinearity leads to a much slower decay (in propagation distance) of the heavy-tailed intensity distribution towards exponential behavior. For stronger nonlinearities even a second, much wider, scintillation peak appears to form. Its height can for very strong nonlinearities even surpass the caustic peak (see figure 13). However, for small disorder strength of the diffracting random potential it became very clear that there is an ideal balance of disorder and nonlinearity for the creation of rogue waves. Whether this is true for stronger disorder we can not conclude from our numerical results so far, but we see indications. More and Figure 11. Probability density of the intensity for different g in potentials with ò=0.0014 (left) and ò=0.014 (right) at a distance L=x c (top panel) and L x 10 c = . As a reference an exponential distribution is shown in purple. Self-focusing nonlinearities enhance the tails of the pdf, while defocusing nonlinearities suppress them. longer simulations for even higher nonlinearities would be needed to answer this questions. If for such strong nonlinearities the paraxial approximation is still reasonable is another question future work will have to answer. Our main conclusion, however, is very clear: to understand the statistics of rogue waves in the presence of disorder (which for example is always there in the ocean in the form of current fields and wind fluctuations) one imperatively needs to understand the interplay of random and nonlinear focusing in detail. We are only interested in y = 0 since the first caustic is on the central line, due to the symmetry of the wavefront and the potential. To obtain an analytically solvable equation we make the approximation that the curvature A t y , 0 0 is constant. Note that ñ is an unspecified constant and depends implicitly on the initial wavefront, i.e. the parameters r, α and β. In the manuscript it is left as a fitting parameter. For g<0 the stability matrix is given by Figure 13. Scintillation index for different ò for strong nonlinearities. For stronger g the increase after the first peak in the scintillation index is also present for larger ò. we get equation (6). In the same way x p g , see equation (7), is derived from 0 This derivation gives a simple functional form of x c g , but does not give quantitative predictions. This is a result of the approximation A t y , 0 const.