Loss of phase and universality of stochastic interactions between laser beams

TL;DR

This paper demonstrates that nonlinear propagation causes laser beams to lose initial phase information, leading to universal interaction statistics that can be modeled without precise phase knowledge.

Contribution

It reveals the loss of phase and universality in laser beam interactions due to nonlinear propagation, introducing a polynomial-chaos approach for statistical modeling.

Findings

Phase information is lost over long distances in nonlinear propagation.

Interaction outcomes become statistically universal, independent of initial phases.

Polynomial-chaos method efficiently predicts interaction statistics.

Abstract

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.

Equations

$$i\frac{\partial }{\partial z}\psi (z,\mathbf{x})+{\nabla }^2\psi +N(|\psi |)\psi =0,$$

$$i\frac{\partial }{\partial z}\psi (z,\mathbf{x})+{\nabla }^2\psi +|\psi {|}^2\psi -ϵ|\psi {|}^4\psi =0,$$

$$i\frac{\partial }{\partial z}\psi (z,\mathbf{x})+{\nabla }^2\psi +\frac{|\psi {|}^2}{1+ϵ|\psi {|}^2}\psi =0.$$

$$\psi (z=0,\mathbf{x};\alpha )={\psi }_0(\mathbf{x};\alpha ),$$

$$\tilde{\phi }(z;\alpha ):=\phi (z;\alpha )\mathrm{mod}(2\pi ).$$

$${\psi }_0(x;\alpha )=3.4(1+0.1\alpha )e^{-x^2},\alpha \sim U(-1,1),$$

$$\psi (z,\mathbf{x};\alpha )\approx e^{i{\eta }_0(\alpha )}{e}^{i\kappa z}{R}_{\kappa }(\mathbf{x})+\mathrm{radiation},$$

$$\phi (z;\alpha )\approx {\eta }_0(\alpha )+z\kappa (\alpha ).$$

$$\Delta \phi (z):=\phi (z;{\alpha }_{\max })-\phi (z;{\alpha }_{\min })\approx \Delta {\eta }_0+z\cdot \Delta \kappa ,$$

$$Z_{\mathrm{lop}}:=\frac{2\pi }{\Delta \kappa }.$$

$${\psi }_0(x)=e^{i\theta x}{R}_{{\kappa }_0}(x-a)+e^{i{\eta }_0}{e}^{-i\theta x}{R}_{{\kappa }_1}(x+a),$$

$$\psi (z,x)\approx e^{i{\kappa }_{0}z}{e}^{i\theta x-i{\theta }^{2}z}{R}_{{\kappa }_0}(x-a-2\theta z)$$

$$+e^{i{\eta }_0}{e}^{i{\kappa }_{1}z}{e}^{-i\theta x-i{\theta }^{2}z}{R}_{{\kappa }_1}(x+a+2\theta z).$$

$$\Delta \phi \approx ({\kappa }_1-{\kappa }_0)z_{\mathrm{cross}}+{\eta }_0.$$

$${\psi }_0(x)=e^{i\theta x}{R}_{{\kappa }_0}(x-a)+c\cdot e^{i{\eta }_0}{e}^{-i\theta x}{R}_{{\kappa }_1}(x+a)$$

$$\text{mbox}\mathbf{random}\kappa :{\eta }_0=0,c=1,{\kappa }_1={\kappa }_0\left(1+\frac{\alpha }{8}\right),$$

$$\text{mbox}\mathbf{random}\kappa \text{mbox},\mathbf{outofphase}:{\eta }_0=\pi ,c=1,$$

$${\kappa }_1={\kappa }_0\left(1+\frac{\alpha }{8}\right),$$

$$\text{mbox}\mathbf{randompower}:{\eta }_0=0,c=1+0.1\alpha ,{\kappa }_1={\kappa }_0,$$

$$\text{mbox}\mathbf{randomphase}:{\eta }_0=\pi \alpha ,c=1,{\kappa }_1={\kappa }_0,$$

$$\underset{z\to \infty }{\lim }\mu \left({\tilde{\phi }}^{-1}\left(\left[a,b\right]\right)\right)=\frac{b-a}{2\pi },$$

$${\phi }_z(x_k^z)=\frac{2k\pi +a}{z},{\phi }_z(y_k^z)=\frac{2k\pi +b}{z},k\in \mathbb{Z}.$$

$$\mu \left({\tilde{\phi }}^{-1}\left(\left[a,b\right]\right)\right)=\sum _{k=k_{\min }}^{k_{\max }}\mu \left(x_k^z,y_k^z\right)+E(z)=$$

$$=\sum _{k=k_{\min }}^{k_{\max }}\mu \left({\phi }_z^{-1}\left(\frac{2\pi k+a}{z}\right),{\phi }_z^{-1}\left(\frac{2\pi k+b}{z}\right)\right)+E(z),$$

$$E(z)=\mu \left([{\alpha }_{\min },y_{k_{\min }-1}^z]\right)+\mu \left([x_{k_{\max }+1}^z,{\alpha }_{\max }]\right).$$

$$z{\phi }_z({\alpha }_{\min }+\delta )-z{\phi }_z({\alpha }_{\min })=$$

$$\left({\eta }_0({\alpha }_{\min }+\delta )-{\eta }_0({\alpha }_{\min })\right)+z\left(\kappa ({\alpha }_{\min }+\delta )-{\kappa }_{\min }\right)$$

$${\alpha }_{\min }\le \underset{z\to \infty }{\lim }{x}_{k_{\min }}^z<{\alpha }_{\min }+\delta .$$

$$\mu \left({\phi }_z^{-1}\left(\frac{2\pi k+a}{z}\right),{\phi }_z^{-1}\left(\frac{2\pi k+b}{z}\right)\right)=$$

$$\underset{{\phi }_z^{-1}\left(\frac{2\pi k+a}{z}\right)}{\overset{{\phi }_z^{-1}\left(\frac{2\pi k+b}{z}\right)}{\int }}c(\alpha )d\alpha =\underset{\frac{2\pi k+a}{z}}{\overset{\frac{2\pi k+b}{z}}{\int }}{g}_z(y)dy,$$

$$\mu \left({\phi }_z^{-1}\left(\frac{2\pi k+a}{z}\right),{\phi }_z^{-1}\left(\frac{2\pi k+b}{z}\right)\right)=$$

$$g_z\left(\frac{2\pi k+{\xi }_k^z}{z}\right)\frac{b-a}{z}.$$

Full text

Nonlinear interactions between two or more laser beams, pulses, and filaments segev1999optical are related to applications ranging from modulation methods in optical communication skidin2016mitigation, to coherent combination of beams fan2005laser; daniault2015xcan; mourou2014fiber; mourou2013future; bellanger2010collective, interactions between filaments in atmospheric propagation couairon2007femtosecond and ignition of nuclear fusion using up to 192 beams kruer1996energy. 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 zakharov1973interaction; ablowitz2004NLS; mitschke1987experimental. In the non-integrable case, however, richer dynamics are possible, including beam repulsion, breakup, fusion and spiraling segev1999optical; snyder1993collisions; Break2001Fusion; tikhonenko1996three. Since the outcome of the interaction strongly depends on the relative phases of the beams garcia1997observation, one can use the initial phase to control the interaction dynamics ishaaya2007self. Nonlinear interactions between solitary waves were also studied in other physical systems scott1973soliton; kivshar1989dynamics, such as fiber optics agrawal2007nonlinear; gordon1983interaction, waveguide arrays meier2004discrete, water waves craik1988wave; su1980head, plasma waves zabusky1965interaction and Bose-Einstein condensates nguyen2014collisions.

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 fibich2011continuations; Merle-92. Therefore, if two such beams intersect after they collapsed, one cannot predict whether they will intersect in- or out-of-phase, and so post-collapse interactions between beams become ”chaotic” and cannot be controlled gaeta2012loss. Loss of phase can also interfere with imaging in nonlinear medium goy2013imaging; barsi2009imaging. 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\pi]$ 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\pi]$, 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.

The propagation of laser beams in a homogeneous medium is governed by the dimensionless nonlinear Schrödinger equation (NLS) in $d+1$ dimensions

[TABLE]

where $z$ is the propagation distance, ${\bf x}=(x_{1},\ldots,x_{d})$ are the transverse coordinates (and/or time in the anomalous regime), and $\nabla^{2}=\partial^{2}_{x_{1}}+\cdots+\partial^{2}_{x_{d}}$ fibich2015nonlinear. Here we consider nonlinearities that support stable solitary waves $\psi=e^{i\kappa z}R_{\kappa}({\bf x})$, such as the cubic-quintic NLS

[TABLE]

or the saturated NLS

[TABLE]

In a physical system the input beam changes from shot to shot. To model this, we write

[TABLE]

where $\alpha$ is the noise realization. We denote by $\varphi(z;\alpha):={\rm arg}\,\psi(z,x=~{}0;\alpha)$ the cumulative on-axis phase at $z$, and study the evolution (in $z$) of the probability distribution function (PDF) of the non-cumulative on-axis phase

[TABLE]

For example, consider the one-dimensional cubic-quintic NLS (2) with the Gaussian initial condition with a random power

[TABLE]

where $U(-1,1)$ is the uniform distribution in $(-1,1)$. Here we consider the one-dimensional case to emphasize that loss of phase and ”chaotic” interactions are not limited to collapsing beams, as was implied by earlier studies fibich2011continuations. At $z=0.15$, the maximal variation of the cumulative phase $\Delta\varphi:=\varphi(\alpha=1)-\varphi(\alpha=-1)$ is fairly small ($\approx~{}0.08\pi$), see Fig. 1($a_{1}$). The corresponding non-cumulative on-axis phase $\tilde{\varphi}$ is identical (Fig. 1($a_{2}$)), and so the probability distribution function (PDF) of $\tilde{\varphi}$, denoted by $f(\tilde{\varphi})$, is fairly localized (Fig. 1($a_{3}$)). As the beam continues to propagate ($z=3$), the maximal variation of the cumulative phase increases to $\Delta\varphi\approx 1.8\pi$, and so $\tilde{\varphi}$ attains most values in $\left[0,2\pi\right]$, though not with the same probability (Fig. 1($b_{1}$)–($b_{3}$)). At an even larger propagation distance ($z=11$), the maximal phase variation is $\Delta\varphi\approx 6.5\pi$, i.e., slightly over three cycles of $\tilde{\varphi}$, see Fig. 1($c_{1}$)–($c_{2}$). At this stage $\tilde{\varphi}$ is nearly uniformly distributed in $[0,2\pi]$, see Fig. 1($c_{3}$), which implies that the beam ”lost” its initial phase $\varphi(z=0)$. By ”loss of phase” we mean that one cannot infer from $f$, the PDF of $\tilde{\varphi}(z;\alpha)$, or from several realizations $\left\{\varphi(z;\alpha_{j}\right\}_{j=1}^{J}$, whether the initial condition was $\psi_{0}(x)=c(\alpha)e^{-x^{2}}$ or $\psi_{0}(x)=e^{i\theta}c(\alpha)e^{-x^{2}}$ for some $0<\theta<2\pi$. Loss of phase is not accompanied by a ”loss of amplitude”. Indeed, the differences between the profiles for $\alpha=\pm 1$ remain small throughout the propagation (Fig. 1($a_{4}$)–($c_{4}$)).

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 $\tilde{\varphi}$ converge to the uniform distribution $f_{\rm uni}(y):\,\equiv\frac{1}{2\pi}$ on $[0,2\pi]$, we plot in Fig. 3(a) the distance $\|f-f_{\rm uni}\|:\,=\int\limits_{0}^{2\pi}\left|f(y)-\frac{1}{2\pi}\right|\,dy$. The convergence is much faster for $d=2$ than for $d=1$, for reasons that will be clarified later 111The convergence is not monotone, because the distance has a local minimum whenever $\Delta\varphi=2\pi k$ for an integer $k$..

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 so

[TABLE]

where $\eta_{0}(\alpha)$ is the on-axis phase which is accumulated during the initial transient, $\kappa$ is the propagation constant, and $R_{\kappa}$ is the positive solution of $\nabla^{2}R_{\kappa}-\kappa R_{\kappa}({\bf x})+N(R_{\kappa})R_{\kappa}=0$. By (7),

[TABLE]

Thus, the nonlinear phase shift grows linearly with $z$ at the rate of $\kappa=\kappa(\alpha)$, see figure 4(b).

Since $\alpha$ is randomly distributed, then so is $\kappa(\alpha)$. More generally, for any initial noise, the beam core evolves into a solitary wave with a random propagation constant $\kappa(\alpha)$, and so $\varphi$ is given by (8) 222This is also the case with multi-parameter noises, e.g., when the amplitude, the phase and the tilt angle are all random.. Consequently, the initial on-axis phase is completely lost as $z\to\infty$:

Proof: see SM.

Lemma 1 provides a new road to the emergence of loss of phase. Indeed, in previous studies fibich2011continuations; gaeta2012loss; Merle-92, 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 $\eta_{0}$ in $\alpha$). Briefly, when a beam undergoes collapse, then in the absence of a collapse-arresting mechanism, $\varphi_{0}(\alpha)\to\infty$ as $z\to Z_{c}(\alpha)$, where $Z_{c}$ is the collapse point fibich2011continuations. Therefore, if a beam undergoes a considerable self-focusing before its collapse is arrested, then it accumulates significant SPM, i.e., $\eta_{0}(\alpha)\gg 2\pi$. In that case, although small changes in $\alpha$ lead to small relative changes in $\eta_{0}(\alpha)$, those are $O(1)$ absolute changes in $\eta_{0}(\alpha)$. In this study, however, we consider non-collapsing beams of the 1D NLS, or mildly-collapsing beams of the 2D NLS. Therefore $\Delta\eta_{0}:=~{}\eta_{0}(\alpha_{\max})-~{}\eta_{0}(\alpha_{\min})\ll~{}2\pi$. 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.

The loss of phase is a nonlinear phenomenon. Indeed, in the linear propagation regime, $\psi(z,{\bf x})=(2iz)^{-\frac{1}{2}}e^{u\frac{|{\bf x}|^{2}}{4z}}\ast\psi_{0}({\bf x})$. Therefore if $\psi_{0}(\alpha_{1})-\psi_{0}(\alpha_{2})\ll 1$ then $\psi(\alpha_{1})-\psi(\alpha_{2})\ll 1$ as well.

Lemma 1 is reminiscent of classical results in ergodic theory of irrational rotations of the circle walters2000ergodic. 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 $\kappa$.

By (8), the maximal difference in the cumulative phase between solutions grows linearly in $z$, i.e.,

[TABLE]

where $\Delta\kappa:\,=\kappa(\alpha_{\max})-\kappa(\alpha_{\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, $\tilde{\varphi}$ is close to be uniformly distributed once $\Delta\varphi(z)\gg~{}2\pi$. Therefore, the characteristic distance for loss of phase is

[TABLE]

Typically, $\Delta\kappa$ is much larger in 2d than in 1d. For example, in Fig. 3(b) $\Delta\kappa\approx 25$ in 2d, and $\Delta\kappa\approx 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 333Denote the solitary-wave power by $P(\kappa):=~{}\int|R_{\kappa}|^{2}\,d{\bf x}$. When $\epsilon=0$ in (2) or (3), then $\frac{dP}{d\kappa}=0$ for $d=2$ but $\frac{dP}{d\kappa}>0$ for $d=1$fibich2015nonlinear. Hence, if $0<\epsilon\ll 1$, then $\frac{dP}{d\kappa}=O(\epsilon)$ for $d=2$, but $\frac{dP}{d\kappa}=O(1)$ for $d=1$. Therefore $\frac{d\kappa}{dP}=O(1)$ for $d=1$ but $\frac{d\kappa}{dP}=O(\frac{1}{\epsilon})$ for $d=2$.. 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 $\psi\to e^{i\beta}\psi$. When the NLS (1) is non-integrable, however, the relative phase between two beams gaeta2012loss; garcia1997observation; krolikowski1997fusion; tikhonenko1996three or condensates nguyen2014collisions 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

[TABLE]

where $a=12$, $\theta=\frac{\pi}{8}$, $\kappa_{0}=8$, and $R_{\kappa}$ is the solitary wave of (2). By Galilean invariance, before the beams intersect at $(z_{\rm cross},x_{\rm cross})\approx(14.7,0)$, each beam propagates as a solitary wave, and so

[TABLE]

[TABLE]

Hence, the difference between the on-axis phases of the two beams at $(z_{\rm cross},x_{\rm cross})$ is

[TABLE]

When the two input beams are in-phase ($\eta_{0}=0$) and identical ($\kappa_{0}=\kappa_{1}$), they intersect in-phase ($\Delta\varphi=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 ($\kappa_{1}=\kappa_{0}+0.1$), then by (11), $\Delta\varphi\approx~{}0.1\cdot 14.7\approx~{}0.48\pi$. 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 $\eta_{0}\approx-0.48\pi$, so that $\Delta\varphi\approx 0$ at $(z_{\rm cross},x_{\rm cross})$, and indeed observe that the two beams merge, see Fig. 5(c) 444Unlike 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..

In what follows, we consider interactions between the two crossing beams

[TABLE]

with four different noise sources:

[TABLE]

where $\alpha\sim U(-1,1)$. Fig. 6(${\text{a}_{1}}$)–(${\text{d}_{1}}$) shows the ”exit intensity” $|\psi(z_{f},x;\alpha)|^{2}$ at $z_{f}=17$ as a function of $x$, for $-1\leq\alpha\leq 1$. As in Fig. 5, depending on $\alpha$, 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_{\rm cross},x_{\rm 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 $\kappa(\alpha)$ is a random variable. Therefore, by Lemma 1, for $z_{\rm cross}\gg Z_{\rm lop}$, the phase of each beam at $(z_{\rm cross},x_{\rm cross})$, and hence also the relative phase between them, is uniformly distributed in $[0,2\pi]$. 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\%\pm 2\%$ and for two filaments is $78\%\pm 2\%$, see Fig. 6(${\text{a}_{2}}$)–(${\text{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\pi]$, 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/\sqrt{N})$ accuracy, where $N$ is the number of NLS solutions. To efficiently use the universal model, we developed a polynomial-chaos 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 phase-dependent methods in optical communications such as Quadrature Amplitude Modulation (QAM) skidin2016mitigation, and in coherent combining of hundreds of laser beams in a small space for ignition of nuclear fusion (Break2001Fusion), and for creating a more powerful laser beam mourou2013future. 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 gauri2017polarization.

Appendix A Proof of Lemma 1

For a given $z\geq 0$, denote $\varphi_{z}(\alpha)=\kappa(\alpha)+\frac{\eta_{0}(\alpha)}{z}$, then $\tilde{\varphi}(\alpha)=z\varphi_{z}(\alpha)\,{\rm mod}\,(2\pi)$. We prove that $\lim\limits_{z\to\infty}\tilde{\varphi}(\alpha)\sim U([0,2\pi])$ by showing that for every $0\leq a<~{}b\leq~{}2\pi$,

[TABLE]

where $\tilde{\varphi}\left[a,b\right]:\,=\left\{\alpha\in[\alpha_{\min},\alpha_{\max}]\,|\,\tilde{u}_{z}(\alpha)\in[a,b]\right\}$.

We first prove the lemma for a strictly monotone function $\kappa$ on $(\alpha_{\min},\alpha_{\max})$. For sufficiently large $z$, $\varphi_{z}$ is also monotone. Let $x_{k}^{z}$ and $y_{k}^{z}$ be the solutions of

[TABLE]

There exists $k_{\min}\leq k_{\max}$ such that $x^{z}_{k_{\min}-1}$ and $y^{z}_{k_{\max}+1}$ do not exist, and for clarity we suppressed the dependence of $k_{\min}$ and $k_{\max}$ on $z$. By definition,

[TABLE]

[TABLE]

where the error term $E(z)$ exists if either $y_{k_{\min}-1}^{z}>\alpha_{\min}$ or $x_{k_{\max}+1}^{z}<\alpha_{\max}$ exist. In such cases, since $d\mu$ is continuous,

[TABLE]

We now show that if $y_{k_{\min}-1}^{z}$ exists, then $\lim\limits_{z\to\infty}\mu\left(\alpha_{\min},y_{k_{\min}-1}^{z}\right)=0$ (a similar proof holds also for $x_{k_{\max}+1}^{z}$). It is enough to show that $\lim\limits_{z\to\infty}x_{k_{\min}}^{z}=\alpha_{\min}$, because $y_{k_{\min-1}}^{z}<x_{k_{\min}}$ and $\mu$ is a continuous measure. Let $\delta>0$, then

[TABLE]

[TABLE]

goes to infinity as $z\to\infty$. Therefore, for large enough $z$, $x_{k_{\min}}^{z}\in(\alpha_{\min},\alpha_{\min}+~{}\delta)$. Thus, for all $\delta>0$,

[TABLE]

Since $\varphi_{z}$ is strictly monotone, then by the inverse function theorem $\varphi_{z}^{-1}\in~{}C^{1}$, and so by substituting $\alpha=\varphi_{z}^{-1}(y)$,

[TABLE]

[TABLE]

where $g_{z}(y):\,=c(\varphi_{z}^{-1}(y))(\varphi_{z}^{-1})^{\prime}(y)$. By Lagrange mean-value theorem, for each index $k$, there exists $\xi_{k}^{z}\in(a,b)$ such that

[TABLE]

[TABLE]

Substituting the above into (15) yields

[TABLE]

Next, consider the integrals

[TABLE]

Denoting $\varphi_{z}(\alpha_{\min}):\,=\varphi_{z,\min}$ and $\varphi_{z}(\alpha_{\max}):\,=\varphi_{z,\max}$ Since

[TABLE]

then $I=I_{2}+O(z^{-1})$. Equating (17b) and (17a), and substituting into (16), yields

[TABLE]

by which we prove (13)

Finally, if $\kappa$, hence $\varphi_{z}$ is piece-wise monotone, we apply the above proof for each sub-interval over which $\varphi_{z}$ is monotone, and by additivity of measure have the result.

Appendix B 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 $\psi(z,x;\alpha)$ 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(\psi)$, e.g., the average intensity over many shots $\mathbb{E}_{\alpha}[\left|\psi\right|^{2}]$.

The standard numerical method for this problem is Monte-Carlo, in which one draws $N$ random values of $\alpha$ and approximates $E_{\alpha}\left[g(\alpha)\right]\approx~{}\frac{1}{N}\sum\limits_{n=1}^{N}g(\alpha_{n})$. The main drawback of this method is its slow $O(1/\sqrt{N})$ convergence rate, where $N$ is the number of NLS simulations. If $g(\alpha):\,=g(\psi(\cdot;\alpha))$ is smooth in $\alpha$, however, we can use orthogonal polynomials as a spectrally accurate basis for interpolation canuto1982approx and numerical integration. Let $\alpha$ is distributed in $[\alpha_{\min},\alpha_{\max}]$ according to a PDF $c(\alpha)$, and let $\left\{p_{n}(\alpha)\right\}_{n=0}^{\infty}$ be the corresponding sequence of orthogonal polynomials, in the sense that $\int\limits_{\alpha_{\min}}^{\alpha_{\max}}p_{n}(\alpha)p_{m}(\alpha)c(\alpha)\,d\alpha=\delta_{n,m}$. For example, if $\alpha$ is uniformly distributed in $[-1,1]$, then $\{p_{n}\}$ are the Legendre polynomials, and if $\alpha$ is normally distributed in $(-\infty,\infty)$, then $\{p_{n}\}$ are the Hermite polynomials. Recall that for smooth solutions one has the spectrally accurate quadrature formula $E_{\alpha}\left[g(\alpha)\right]\approx\sum\limits_{j=1}^{N}g(\alpha_{j}^{N})w_{j}^{N}\,,$ where $\{\alpha_{j}^{N}\}_{j=1}^{N}$ and $\left\{w_{j}^{N}\right\}_{j=1}^{N}$ are the roots of the orthogonal polynomial $p_{N}(\alpha)$ and their respective weights $w_{j}^{N}=\int\limits_{\alpha_{\min}}^{\alpha_{\max}}\prod\limits_{i=1,\,i\neq j}^{N}\frac{\alpha-\alpha_{i}^{N}}{\alpha_{j}^{N}-\alpha_{i}^{N}}\,c(\alpha)\,d\alpha$ 555See day2005roots for a numerically efficient and stable algorithm for computing the $\left\{\alpha_{j}^{N},\,w_{j}^{N}\right\}_{j=1}^{N}$.. We apply the collocation Polynomial Chaos Expansion (PCE) method as follows o2013polynomial; xiu2010numerical:

For $j=1,\dots,N$, solve the NLS for $\psi\left(z,{\bf x};\alpha_{j}^{N}\right)$, and set $g(\alpha_{j}^{N}):=g\left(\psi(z,x;\alpha_{j}^{N})\right)$. 2. 2.

Approximate

[TABLE]

This method is ”non-intrusive”, i.e., it does not require any changes to the deterministic NLS solver. Moreover, the orthogonality of $\{p_{n}\}$ leads to direct formulae for the mean and standard deviation of $g$:

[TABLE]

[TABLE]

As noted, the PCE method has a spectral convergence rate for smooth functions. For example, the results in Fig. 1 2 were computed using $N=10$ and $N=31$ NLS simulations, respectively. To reach a similar accuracy with the Monte Carlo method would require more than $1000$ NLS simulations. Some quantities of interest, however, such as the number of filaments (Fig. 6($a_{2}$)–6($d_{2}$)), or the non-cumulative on-axis phase $\tilde{\varphi}={\rm arg}\,\left(\psi(z,{\bf x}=0;\alpha)\right)\,{\rm mod}(2\pi)$, (Fig. 1($a_{3}$)–1($c_{3}$) and Fig. 2($a_{3}$)–2($c_{3}$)) are non-smooth. Therefore, a straightforward application of the PCE method for such quantities requires $O\left(10^{3}\right)$ simulations to converge. In such cases, we begin with stages (1)–(2) and calculate the PCE approximation (18) of the smooth function $\psi(z,{\bf x};\alpha)$ using $\{\psi(z,{\bf x};\alpha_{j}^{N})\}_{j=1}^{N}$ with a relatively small $N$. Then we proceed as follows:

Use the gPC interpolant (18) to obtain $\psi(\cdot,\tilde{\alpha}_{m})\approx\psi_{N}(\cdot,\tilde{\alpha}_{m})$ on a sufficiently dense grid $\{\tilde{\alpha}_{m}\}_{m=1}^{M}$, where $M\gg N$. 2. 4.

Compute $g(\tilde{\alpha}_{m})\approx g(\psi_{N}(\cdot,\tilde{\alpha}_{m}))$ for $m=1,\ldots,M$, 3. 5.

Compute the statistics of $g(\psi)$ using $\left\{g(\psi_{N}(\cdot,\tilde{\alpha}_{m}))\right\}_{m=1}^{M}$.

For example, when we computed the number of beams at $z=z_{f}$ in Fig. 6, we first computed the PCE interpolant $\psi_{N}(z_{f},x;\alpha)$ with $N=71$. Then we computed $\psi(z_{f},x;\tilde{\alpha}_{m})\approx\psi_{N}(z_{f},x;\tilde{\alpha}_{m})$ for $m=1,\ldots,M=801$. For each $\tilde{\alpha}_{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 $\psi_{N}$ (18) at $M\gg N$ grid points in step (3) is negligible compared to directly solving the NLS for $N$ times in stage (1).