Loss of phase and universality of stochastic interactions between laser beams

We show that all laser beams gradually lose their initial phase information in nonlinear propagation. Therefore, if two beams travel a sufficiently long distance before interacting, it is not possible to predict whether they would intersect in- or out-of-phase. Hence, if the underlying propagation model is non-integrable, deterministic predictions and control of the interaction outcome become impossible. Because the relative phase between the two beams becomes uniformly distributed in $[0,2\pi]$, however, the statistics of the interaction outcome are universal, and can be efficiently computed using a polynomial-chaos approach, even when the distributions of the noise sources are unknown.

We show that all laser beams gradually lose their initial phase information in nonlinear propagation. Therefore, if two beams travel a sufficiently long distance before interacting, it is not possible to predict whether they would intersect in-or out-of-phase. Hence, if the underlying propagation model is non-integrable, deterministic predictions and control of the interaction outcome become impossible. Because the relative phase between the two beams becomes uniformly distributed in [0, 2π], however, the statistics of the interaction outcome are universal, and can be efficiently computed using a polynomial-chaos approach, even when the distributions of the noise sources are unknown. Nonlinear interactions between two or more laser beams, pulses, and filaments [1] are related to applications ranging from modulation methods in optical communication [2], to coherent combination of beams [3][4][5][6][7], interactions between filaments in atmospheric propagation [8] and ignition of nuclear fusion using up to 192 beams [9]. In the integrable one-dimensional cubic case, such interactions can only lead to phase and lateral shifts, which can be computed analytically using the Inverse Scattering Transform [10][11][12]. In the non-integrable case, however, richer dynamics are possible, including beam repulsion, breakup, fusion and spiraling [1,[13][14][15]. Since the outcome of the interaction strongly depends on the relative phases of the beams [16], one can use the initial phase to control the interaction dynamics [17]. Nonlinear interactions between solitary waves were also studied in other physical systems [18,19], such as fiber optics [20,21], waveguide arrays [22], water waves [23,24], plasma waves [25] and Bose-Einstein condensates [26].
In previous studies it was shown, both theoretically and experimentally, that when a laser beam undergoes an optical collapse, its initial phase is "lost", in the sense that the small shot-to-shot variations in the input beam lead to large changes in the nonlinear phase shift of the collapsing beam [27,28]. Therefore, if two such beams intersect after they collapsed, one cannot predict whether they will intersect in-or out-of-phase, and so postcollapse interactions between beams become "chaotic" and cannot be controlled [29]. Loss of phase can also interfere with imaging in nonlinear medium [30,31]. Note that loss of phase does not imply a loss of coherence, but rather that at any given propagation distance, the coherent beam can only be determined up to an unknown constant phase.
In this study we show that loss of phase is ubiquitous in nonlinear optics. Thus, while collapse accelerates the loss of phase process, non-collapsing or mildly-collapsing beams can also undergo a loss of phase. The loss of phase builds up gradually with the propagation distance, i.e., the shot-to-shot variations of the beam's nonlinear phase shift increase with the propagation distance, and approach a uniform distribution in [0, 2π] at sufficiently large distances. As noted, because of the loss of phase, deterministic predictions and control of interactions between laser beams become impossible. We show, however, that loss of phase allows for accurate predictions of the statistical properties of these stochastic interactions, even without any knowledge of the noise source and characteristics. Indeed, because the relative phase between the beams becomes uniformly distributed in [0, 2π], the statistics of the interaction are universal, and can be computed using a "universal model" in which the only noise source is a uniformly distributed phase difference between the input beams. These computations can be efficiently performed using a polynomial-chaos based approach.
We obtained similar results for the same equation and initial condition in two dimensions, see Fig. 2. To compare the rates at which the PDFs ofφ converge to the uniform distribution f uni (y) : ≡ 1 2π on [0, 2π], we plot in Fig. 3 The convergence is much faster for d = 2 than for d = 1, for reasons that will be clarified later [40].
To understand the emergence of loss of phase in Fig. 1-2, we note that after an initial transient, the beam core evolves into a stable solitary wave, see e.g., Fig. 4(a), and where η 0 (α) is the on-axis phase which is accumulated during the initial transient, κ is the propagation constant, and R κ is the positive solution of (7), Thus, the nonlinear phase shift grows linearly with z at the rate of κ = κ(α), see figure 4(b).
Since α is randomly distributed, then so is κ(α). More generally, for any initial noise, the beam core evolves into a solitary wave with a random propagation constant κ(α), and so ϕ is given by (8) [41]. Consequently, the initial on-axis phase is completely lost as z → ∞: Lemma 1. Let α be a random variable which is distributed in [α min , α max ] with an absolutelycontinuous measure dµ, let κ(α) be a continuously differentiable, piece-wise monotone function on [α min , α max ], let η 0 (α) be continuously differentiable on [α min , α max ], and let ϕ be given by (8).
Proof: see SM. Lemma 1 provides a new road to the emergence of loss of phase. Indeed, in previous studies [27][28][29], the loss of phase was caused by the large self-phase modulations (SPM) that accumulate during the initial beam collapse (i.e., by the variation of η 0 in α). Briefly, when a beam undergoes collapse, then in the absence of a collapsearresting mechanism, ϕ 0 (α) → ∞ as z → Z c (α), where Z c is the collapse point [27]. Therefore, if a beam undergoes a considerable self-focusing before its collapse is arrested, then it accumulates significant SPM, i.e., η 0 (α) ≫ 2π. In that case, although small changes in α lead to small relative changes in η 0 (α), those are O(1) absolute changes in η 0 (α). In this study, however, we consider non-collapsing beams of the 1D NLS, or mildly-collapsing beams of the 2D NLS. Therefore ∆η 0 := η 0 (α max ) − η 0 (α min ) ≪ 2π. In such cases, the loss of phase builds up gradually with the propagation distance z, and not abruptly during the initial collapse, as in the previous studies.
Lemma 1 is reminiscent of classical results in ergodic theory of irrational rotations of the circle [33]. Unlike these results, however, Lemma 1 does not describe the trajectory of a single point on the circle under consecutive discrete phase additions, but rather the convergence of a continuum of points under with continuous linear change with a varying rate κ.
By (8), the maximal difference in the cumulative phase between solutions grows linearly in z, i.e., where ∆κ : = κ(α max )−κ(α min ) is the maximal variation in the propagation constant, induced by the noise. As suggested by the proof of Lemma 1 and by Figs. 1 and 2, ϕ is close to be uniformly distributed once ∆ϕ(z) ≫ 2π. Therefore, the characteristic distance for loss of phase is Typically, ∆κ is much larger in 2d than in 1d. For example, in Fig. 3(b) ∆κ ≈ 25 in 2d, and ∆κ ≈ 1.8 in 1d. Intuitively, this is because the input beam evolves into a solitary wave, and over a given power range, the propagation constant of the solitary wave changes considerably less in 1d than in 2d [42]. Hence, by (9), the loss of phase occurs much faster in the two-dimensional case than in the one-dimensional case, thus explaining Fig. 3(a). A priori, the loss of initial phase has no physical implications, since the NLS (1) is invariant under the transformation ψ → e iβ ψ. When the NLS (1) is non-integrable, however, the relative phase between two beams [15,16,29,34] or condensates [26] can have a dramatic effect on their interaction, thus making the loss of initial phase physically important. To illustrate that, consider again the cubic-quintic NLS (2) for d = 1 with the two crossing beams initial condition ψ 0 (x) = e iθx R κ0 (x − a) + e iη0 e −iθx R κ1 (x + a) , (10) where a = 12, θ = π 8 , κ 0 = 8, and R κ is the solitary wave of (2). By Galilean invariance, before the beams intersect at (z cross , x cross ) ≈ (14.7, 0), each beam propagates as a solitary wave, and so Hence, the difference between the on-axis phases of the two beams at (z cross , x cross ) is When the two input beams are in-phase (η 0 = 0) and identical (κ 0 = κ 1 ), they intersect in-phase (∆ϕ = 0), and so they merge into a strong central beam, see Fig. 5(a). If we introduce a 1.25% change in the propagation constant of the second beam (κ 1 = κ 0 + 0.1), then by (11), ∆ϕ ≈ 0.1 · 14.7 ≈ 0.48π. This phase difference is sufficient for the two beams to repel each other, see Fig. 5(b). Therefore, the interaction is "chaotic", in the sense that a small change in the input beams leads to a large change in the interaction pattern. To further demonstrate that the change in the interaction pattern is predominately due to the phase difference, we "correct" the initial phase of the second beam by setting η 0 ≈ −0.48π, so that ∆ϕ ≈ 0 at (z cross , x cross ), and indeed observe that the two beams merge, see In what follows, we consider interactions between the two crossing beams with four different noise sources: where α ∼ U (−1, 1). Fig. 6(a 1 )-(d 1 ) shows the "exit intensity" |ψ(z f , x; α)| 2 at z f = 17 as a function of x, for −1 ≤ α ≤ 1. As in Fig. 5, depending on α, there are two possible outputs: Either a single beam (resulting from beam fusion), or two beams (resulting from beam repulsion). Generally speaking, there is a single output beam whenever the two input beams are "sufficiently" in-phase at (z cross , x cross ). In a physical setting the noise distribution is typically unknown. Nevertheless, the on-axis phase of each beam core is given by (8), where κ(α) is a random variable. Therefore, by Lemma 1, for z cross ≫ Z lop , the phase of each beam at (z cross , x cross ), and hence also the relative phase between them, is uniformly distributed in [0, 2π]. Hence, to leading order, the statistics of the interactions are universal. Indeed, for all 4 noisy initial conditions we observe that: (i) the probability for a single filament is 22% ± 2% and for two filaments is 78% ± 2%, see Fig.  6(a 2 )-(d 2 ), and (ii) the mean and standard deviation of the lateral locations of the output beams are nearly identical, see Fig. 6(e).
The above results show that the statistics of long-range interactions between laser beams are independent of the noise source and its characteristics, and can be computed using a "universal model" in which the only randomness comes from the addition of a random constant phase to one of the input beams, which is uniformly distributed in [0, 2π], c.f. (12d). The standard approach for computing the statistics of the interactions in this "universal model" is the Monte-Carlo method. This method, however, is very inefficient due to its O(1/ √ N ) accuracy, where N is the number of NLS solutions. To efficiently use the universal model, we developed a polynomialchaos based numerical method, which is both spectrally accurate and makes use of any deterministic NLS solver. For further details, see SM.
In conclusion, we showed that when laser beams or pulses interact after a sufficiently long propagation distance, their relative phase at the crossing point cannot be predicted or controlled. In such cases, the notion of a "typical experiment" or a "typical solution" may be misleading, and one should adopt a stochastic approach. The loss of phase can explain some of the difficulties in phasedependent methods in optical communications such as Quadrature Amplitude Modulation (QAM) [2], and in coherent combining of hundreds of laser beams in a small space for ignition of nuclear fusion [14], and for creating a more powerful laser beam [6]. In these applications, controlling the phases of the input beams or pulses might not provide a good control over their interaction or combination, due to the loss of phase. Our study suggests that controlling the relative phases at the intersection point may be achieved by either shortening the propagation distance, or by coupling the beams throughout the propagation. Loss of phase is also relevant to the loss of polarization for elliptically-polarized beams [35]. , then I = I 2 + O(z −1 ). Equating (17b) and (17a), and substituting into (16), yields by which we prove (13) Finally, if κ, hence ϕ z is piece-wise monotone, we apply the above proof for each sub-interval over which ϕ z is monotone, and by additivity of measure have the result.

Numerically solving the universal model
Although in the universal model the noise is uniformly distributed, we allow for a more general noise distribution, so that we can e.g., produce results such as figure 1 for non-uniform noise distributions.
Let ψ(z, x; α) be the solution of the NLS (1) with the random initial condition (6). In what follows, we introduce an efficient numerical method for computing the statistics of g(ψ), e.g., the average intensity over many The standard numerical method for this problem is Monte-Carlo, in which one draws N random values of α and approximates E α [g(α)] ≈ 1 N N n=1 g(α n ). The main drawback of this method is its slow O(1/ √ N ) convergence rate, where N is the number of NLS simulations. If g(α) : = g(ψ(·; α)) is smooth in α, however, we can use orthogonal polynomials as a spectrally accurate basis for interpolation [36] and numerical integration. Let α is distributed in [α min , α max ] according to a PDF c(α), and let {p n (α)} ∞ n=0 be the corresponding sequence of orthogonal polynomials, in the sense that where {α N j } N j=1 and w N j N j=1 are the roots of the orthogonal polynomial p N (α) and their respective weights . We apply the collocation Polynomial Chaos Expansion (PCE) method as follows [37,38]: 1. For j = 1, . . . , N , solve the NLS for ψ z, x; α N j , and set g(α N j ) := g ψ(z, x; α N j ) .
For example, when we computed the number of beams at z = z f in Fig. 6, we first computed the PCE interpolant ψ N (z f , x; α) with N = 71. Then we computed ψ(z f , x;α m ) ≈ ψ N (z f , x;α m ) for m = 1, . . . , M = 801. For eachα m , we count the number of filaments and used this to produce the histogram in figure 6(a 2 )-6(d 2 ). The additional computational cost of sampling ψ N (18) at M ≫ N grid points in step (3) is negligible compared to directly solving the NLS for N times in stage (1).
[41] This is also the case with multi-parameter noises, e.g., when the amplitude, the phase and the tilt angle are all random.
[43] Unlike Fig. 5(a), the output beam is slightly tilted upward, since the lower input beam is more powerful, and therefore the net linear momentum points upward.
[44] See [39] for a numerically efficient and stable algorithm for computing the α N j , w N j N j=1 .