Random walk of solitary and shock waves in nonlocal disordered media

We study the random walk of solitons and characteristic lines of shock fronts in the presence of disorder for the one-dimensional nonlinear Schrödinger equation in Kerr-like media. We analyze the interplay of nonlocality and randomness, and the way their competition affects strongly coherent nonlinear waves is theoretically and numerically investigated.


Introduction
Nonlinear phenomena, as solitons [1,2] and dissipative and dispersive shock waves (SWs) [3][4][5][6][7][8], are largely affected by the presence of nonlocality [9][10][11] and disorder [12][13][14]. A nonlocal response is present when the light beam induces a refractive index change on a spatial region greater than the beam waist; this effect is typically found in physical systems displaying long-range correlation [15][16][17][18][19][20]. However, considering a nonlocal response also implies accounting for disorder. Indeed, a large spatial region that interacts with the propagating electromagnetic field in a nonlocal medium commonly also involves material fluctuations. These are due to electronic and thermal effects and their interplay with nonlinearity produces a complex scenario to be investigated. In fact, the way nonlinear coherent phenomena, as SWs and solitons, are affected by randomness is one of the most active research areas in nonlinear optics [21][22][23][24]. Disorder frustrates the coherent interaction between light and matter and delays the formation of solitons or SWs. In these respects, the amount of disorder and nonlinearity can be retained as control parameters in 'phase-diagrams' [14,25,26], identifying regimes in which various phenomena can be observed.
In optics, one-dimensional (1D) propagation in a nonlinear medium is described by the nonlinear Schrödinger equation (NLS) [1] 2ik∂ z A + ∂ 2 x A + 2k 2 n n 0 where A is the beam envelope normalized such that I = |A| 2 is the optical intensity, n(r ) = n 2 |A| 2 is the intensity-dependent component of the refractive index, n = n 0 + n 2 |A| 2 . z is the propagation distance, x the transverse coordinate, k = 2πn 0 /λ is the wave vector and λ is the wavelength.
In this paper, we introduce a disorder term in (1) by taking n(x, z) = n 2 |A| 2 + n R , where n R = n 0 R V (x, z) and V (x, z) is a Gaussian potential such that n R (x) n R (x ) = ( n 0 R ) 2 δ(x − x ). The refractive index perturbation n R is associated to a length scale competing with the scale of the nonlinear interaction between light and matter, as detailed below.
Equation (1) with the disorder term (2) can be normalized by introducing spatial scales weighting the relative roles of disorder, nonlinearity and diffraction. Nonlocality [9,10,14,23] also introduces a spatial scale. The nonlocal response of the medium can be written as a convolution integral of the wave field and a response function χ(x − x ): The response function typically depends on a length scale L nloc which gives its spatial extension, being the local limit recovered by . The typical response function has an exponential form χ (x) = e −|x|/L nloc /(2L nloc ). In the highly nonlocal limit, the response function can be approximated by a constant, χ (x) = χ 0 , with χ 0 = 1/2L nloc in the exponential case.

The three length scales
When we add disorder and nonlocality in (1), the resulting equation is where we have introduced the convolution integral The diffraction (or Rayleigh) length L d for a Gaussian beam, propagating along the z-axis and with a beam waist w 0 , is defined as The weight of nonlinearity is measured by the nonlinear length L nl , which defines the length over which the nonlinear phase shift φ nl = 4π n 2 I 0 z/λ is equal to 1: with I 0 the peak intensity. We define the randomness length L R as the length over which the phase shift of a beam propagating in an homogeneous medium with refractive index n 0 R , φ R = 4π n 0 R z/λ is equal to 1: and it is determined by the strength of the random refractive index perturbation n 0 R . n 0 R is related to the localization length l, which scales as l ∝ 1/( n 0 R ) 2/3 [12,27], and measures the spatial extension of the disorder induced transversal confinement.

Rescaled nonlinear Schrödinger equation
Scaling the coordinates as x → w 0 x, z → 2Lz with L ≡ √ L d L nl , and by introducing the complex field ψ = A/ √ I 0 , we obtain from (4) where and are the only relevant quantities, evidencing the relative weight of nonlinearity, randomness and diffraction. The rescaled response function in the exponential case takes the form K (x) = 1/(2σ ) e −|x|/σ , with σ = L nloc /w 0 .
In the next section, we show how different regimes can be obtained when varying the length scales L d , L nl and L R , and correspondingly, 1 and 2 .

Ordered case
In an homogeneous sample, the weight of the fluctuations of refractive index n 0 R is negligible, hence the randomness length scale is infinite L R L nl , L d or, equivalently, 2 = 0. The rescaled NLS equation becomes where we consider the local case K (x) = δ(x).

Solitonic regime
The generation of solitary waves results from the balance between the diffraction and nonlinearity, in this case L d Ł nl , 1 1. Limiting ourselves to the focusing case, equation (8) i∂ z ψ + ∂ 2 x ψ + |ψ| 2 ψ = 0 admits a localized solution, taking the form ψ(x, z) = u(x) e iβz with W 0 = √ 2/u 0 and β = u 2 0 /2. The condition on the length scales L d = L nl gives 'the existence curve', u 0 W 0 = const that relates amplitude and waist of soliton. In the 1D case, the localized solution is always stable and no catastrophic collapse events occur. In the nonlocal case, u(x) is found numerically [14].

Shock waves regime
The SWs generation is a strongly nonlinear phenomenon and has also been studied in the presence of nonlocality [9][10][11]. To investigate such a regime, we deal with the hydrodynamical approximation [9,28], when considering a weakly diffractive regime, L d L nl and 1 1. The evolution of a Gaussian beam ψ(x) = u 0 (x) e iβz with u 0 (x) = e −x 2 /2 is obtained by the Wentzel, Kramers and Brillouin (WKB) transformation where ρ(x, z) is the slowly varying beam intensity and φ(x, z) its phase. By inserting equation (14) in (11), the following Euler-like equations are obtained: where we have defined the velocity of the fluid v ≡ φ x as the phase chirp and we neglect the 'quantum' pressure term, 1 2 At the beginning of the evolution, the strong nonlinearity lets the field intensity ρ to act as a pressure distribution on the velocity of the 'optical fluid'. In the focusing case, the pressure profile induces the fluid (the flat phase distribution) to accumulate in the center causing the shock to occur in correspondence of the peak intensity. In the defocusing case, this results in an accumulation of the energy density on the beam sides with the formation of two symmetric shocks, as described below.

Weakly disordered system
When the disorder is a perturbation, L R L nl , L d and ε 2 ε 1 , it is interesting to analyze the way it affects the evolution of solitary and SWs.

Solitonic regime
We consider a nonlocal nonlinear medium in the presence of disorder. We assume that L R L nl , L d and the diffraction and nonlinearity satisfy the condition for soliton existence: L nl L d . As a consequence, equation (8) becomes where we introduced the variable η = ε 2 2 V (x, z) = ( n 0 R /|n 2 |I 0 )V (x, z) as a small perturbation. We use a perturbational approach and show that an increasing degree of nonlocality hampers the Brownian motion of self-trapped wave-packets [14]. In fact, the nonlocality acts as a spatial filter on the wave propagation, canceling fluctuations.
Letting ψ = φ exp(iβz), equation (17) is written as where s is taken as a perturbation term, depending on ψ and its transverse derivative at any order, and β is the nonlinear wave-vector. For s = 0, the unperturbed solitary wave (SW) is written as where X is the center of the self-trapped wave, θ is the phase and 2 is the velocity. In the following, we derive dynamical equations for the first order perturbation of these parameters. By letting φ = φ 0 + φ 1 , the linearized evolution equation is with The analysis can be limited to SW with = 0. The first order perturbation φ 1 is represented by a small variation of X , , θ and β: where the auxiliary functions are defined as f θ = iφ 0 , f β = ∂ β φ 0 , f X = ∂ X φ 0 and f = −i(x − X )φ 0 , and being δ X (z), δθ (z), δβ(z) and δ (z) the z-dependent perturbations to SW parameters. We introduce the adjoint functionsf given byf with a and b two SW parameters (X , , θ or β) and introducing the scalar product BeingL(i f ) = −iL( f ), the projection over the adjoint functions of equation (20) gives where S α = (f α , s), and the dot is the derivative with respect to z. Equations (23) hold for any s; for the random perturbation in (17), we have where from which and averaging over disorder, where From (23), we obtain the mean position from which δ X = 0 and Equation (30) is the nonlocal expression of the Gordon-Haus effect, that describes the Brownian motion of coherently amplified solitons [13]. In equation (30), the quantity C measures the effect on locality on the soliton fluctuations. In particular, in the highly nonlocal limit K (x) → K 0 , for a bell-shaped soliton profile [u(x) = u(−x)], we have In the highly nonlocal regime, the random fluctuations of the fundamental soliton vanish. Physically, when we increase the degree of nonlocality, the spectral content of the nonlocal response shrinks, and, as a consequence, averages out the randomness [15]. This result is valid for every kind of nonlocality.
To verify this theoretical analysis, we have numerically integrated the stochastic partial differential equation (17) for a 1D exponential nonlocality, using a pseudo-spectral stochastic Runge-Kutta algorithm [29,30]. Panels (a) and (b) of figure 1 show a typical evolution starting from a bound state and displaying the random deviation of the SW for two different degrees of nonlocality σ . We stress that the high nonlocality degree for the case in panel (b) considerably dampens the random walk of the soliton with respect to the quasi-local case in panel (a). In panels (c) and (d) of figure 1, we report various trajectories for a fixed SW power for σ 2 = 0.4 (c) and σ 2 = 10 (d). Panel (e) shows the calculated standard deviation for the two degrees of nonlocality above: the analytical prediction of equation (30) is well reproduced by the numerical simulations.

Shock waves regime
Due to their strongly nonlinear and coherent origin, the SWs will be heavily affected by disorder during focusing and defocusing propagation. By the hydrodynamical approach, we investigate the competition between strong nonlinearity and disorder and we find that an increasing amount of disorder delays the shock point up to the inhibition of the wave breaking.
To enter the shock regime, we consider a weak diffraction regime with L d L R L nl , and the local case, K = δ(x): letting ψ(x, z) = √ ρ(x, z) e iφ(x,z)/ 1 , and retaining only leading orders in 1 , we have where v ≡ φ x is the phase chirp. We derive the first order equation with respect to x: where the potential function U is U foc = η − ρ in the focusing case and U def = η + ρ in the defocusing case. We assume a Gaussian profile for the input beam ρ = exp(−x 2 /2), and we apply the method of characteristics [31]. Equation (34) is reduced to dx dz = 2v, and by deriving the first equation with respect to z (hereafter, we make the ansatz z → z/ √ 2 to simplify the notation), we have which is equation of motion of a unitary mass particle in the potential U = U nl + η with U nl = e −x 2 /2 the deterministic contribution deriving from nonlinearity. The stochastic term f = −dη/dx is taken as a Langevin force with Gaussian distribution, Numerical simulations, obtained by a stochastic Runge-Kutta algorithm [32], provide solutions of equation (36) for several initial conditions x(0), v(0) and noise realizations η(x). For each pair of initial conditions, we change the noise configuration in order to account for the dependence of η N on x. The shock position corresponds to the singularity of the phase chirp |dv/dx| → ∞ and can be directly measured as the point of large steepness of the beam intensity [10]. In this unitary mass particle approach and in the absence of disorder, the shock point can be extrapolated by the coalescence of multiple trajectories, as shown in figure 2(a) (defocusing case) and in figure 3(a) (focusing case). In fact, the characteristic lines give the direction of energy flow and the shock formation corresponds to the accumulation of light density in specific regions of the transverse dimension x: at the edges of the laser beam for defocusing nonlinearity, and at the center of the laser beam for focusing nonlinearity. When we analyze the phase space of x and v (figure 2(c) for defocusing nonlinearity and figure 3(c) for focusing nonlinearity), the shock can be observed as a folding of particle velocities into a multivalued function. Physically, being that the light field is a single-value function, it happens that diffraction regularizes the wavefront, breaking the wave and causing the so-called undular bore [8,16,33]. This phenomenon corresponds to coherent fast oscillations both of the field intensity and the phase chirp, which smoothly join the singularity region with the rest of the wave front.
In the presence of disorder, the light scattering lowers the field intensity and moves the shock threshold away along the propagation distance. In the particle analogy, the characteristic lines are randomized (figures 2(b) and 3(b)) and the particles collide along greater distances, undergoing a Brownian motion and delaying the trajectory accumulation that determines the shock event.
In figures 4 and 5, we show, respectively, the numerical histograms of the particle positions in the defocusing and focusing cases. Panel (a) reports the accumulation point of the caustics in the ordered systems; it has a fixed position, corresponding to the singularity of the phase chirp, dv/dx → ∞. In the disordered cases, panel (b), the trajectories are randomized and intersect at different points. We choose as the shock position z s (η N ), the mean value among the positions in which the histograms in figures 4 and 5 (panels (b)) are above the 90% of its maximum. It is an arbitrary criterion, which allows us to assess the probability that the shock event can be found in a specific spatial region.
In panels (c) of figures 4 and 5, the shock position is reported as a function of normalized strength of disorder η N . As η N increases, the shock process is delayed with respect to the ordered case. When the amount of disorder exceeds a critical value, the shock generation is inhibited. This threshold can be analytically determined as that corresponding to a random index perturbation n R comparable with the nonlinear term: ( n 0 R ) 2 1/2 ∼ = |n 2 |I 0 , which gives η N ∼ = 1. In equation (32), the stochastic term related to η becomes greater than the nonlinear contribution. When this happens, the hydrodynamical limit in (36) is no longer valid. In fact, the stochastic fluctuations of the material become so large as to mask the nonlinear effect.
In the hydrodynamical regime, the shock position z s scales as 1/ √ P. When disorder is added, we found that the shock position becomes like z s (η N )/ √ P with z S (η N ) given in figures 4(c) and 5(c). Hence, nonlinearity and disorder act as opposite effects: on the one hand, the shock formation is favored (z S decreases) when nonlinearity grows, but on the other hand, by introducing disorder, light is scattered and the phenomenon is hampered (z s increases).

Nonlocal case.
Here we consider the way nonlocality affects the formation of SWs in the presence of disorder. Similarly to disorder, nonlocality delays the shock generation. Furthermore, nonlocality also allows the shock to form in focusing media whereas filamentation or modulation instability prevail on the shock formation in local media [9].
We start by considering the nonlocal version of equation (32) where θ = 2k L nl n/n 0 = K * (±|ψ| 2 + η) and n is the nonlocal nonlinear disordered refractive index perturbation (3) with a degree of nonlocality given by σ . We make a WKB transformation of the field ψ(x, z) = √ ρ(x, z) e iφ(x,z) 1 and substituting in (37), we obtain for the leading orders in 1 : where v = ∂ x φ is the phase chirp. The intensity ρ adiabatically changes in z with respect to the phase chirp and the refractive index modulation θ . Hence, we take ρ constant versus z, ρ(x) = e −x 2 /2 . As a consequence, the exact solution θ(x) can be obtained from the third equation in (38). By using the method of characteristics, the second equation of (38) can be reduced to dx dz = 2v, where the potential U is U foc = K * (η − ρ) in the focusing case and U def = K * (η + ρ) in the defocusing case. The nonlocal nonlinear response K * ρ is the counterpart of the deterministic potential represented by the field intensity ρ in the local case (see (36)).
In the nonlocal case, the potential well (focusing case) and the potential barrier (defocusing case) are much more flat with respect to the local case. As a consequence, the intersection of characteristic lines leading to the shock generation is slowed down and the shock position is delayed along z. The potential profiles are randomized by the presence of K * η in (39) that further act as a SW damper, corrugating the surface over which the effective particles evolve.
Being the SWs generated in correspondence of the chirp singularity |dφ/dx| → ∞ and, approximating the phase as φ ≈ θ , we have In the highly nonlocal limit K (x − x ) = K 0 , and the integral (40) vanishes. As a result, by increasing the nonlocality degree, the shock position is gradually delayed, up to its inhibition in the highly nonlocal limit. Hence, both nonlocality and disorder tend to destroy the shock formation. The phenomena dampen the effective pressure gradient that causes the thickening of the 'optical fluid' resulting in the shock generation. On the one hand, nonlocality flattens the pressure profile over the spatial region of the nonlocal response (given by σ ). On the other hand, disorder scatters the light field, reducing the field intensity accumulation on the shock points.

Conclusions
We have investigated the interplay between nonlinearity and disorder in the nonlocal NLS. Three main length scales, namely diffraction, randomness and nonlinear lengths, control the excitation of solitary and SWs in weakly disordered media. When solitons are considered, the nonlocality limits the scattering action of disorder on the soliton trajectory, averaging out the randomness and stabilizing the solitary waves. In the SW regime, a nonlocal nonlinearity delays shock formation with respect to the local case. The presence of disorder further delays the shock generation while nonlocality on the random potential does not significantly change the scenario. Given that disorder is inextricably linked to nonlocality, we believe that these results can be relevant when considering thermal and thermodiffusive media [7,9,34], photorefractive systems [17], Bose-Einstein condensate [18,19], liquid crystals [20] or plasma [35,36].