# Effect of atmospheric turbulence on modulational instability in laser-pulse propagation

Alain Bourdier, Jean-Claude Diels, Hassen Ghalila, Olivier Delage

► To cite this version:

Alain Bourdier, Jean-Claude Diels, Hassen Ghalila, Olivier Delage. Effect of atmospheric turbulence on modulational instability in laser-pulse propagation. Physica D: Nonlinear Phenomena, 2024, 457, 10.1016/j.physd.2023.133974. insu-04846365

HAL Id: insu-04846365

https://insu.hal.science/insu-04846365v1

Submitted on 18 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Effect of atmospheric turbulence on modulational instability in laser-pulse propagation

Alain Bourdier<sup>a, *</sup>, Jean-Claude Diels<sup>a</sup>, Hassen Ghalila<sup>b</sup>, Olivier Delage<sup>c</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

Communicated by Prof Feng Bao-Feng

**Keywords:**

Laser pulse propagation

Nonlinear optics

Multiple filamentation

Instability

Turbulence

Stochastic differential equations

Kubo-Anderson process

## ABSTRACT

Multiple filamentation is a major problem for laser pulse propagation in the atmosphere. In this article, we study the influence of a turbulent atmosphere on the growth of the modulational instability which is the cause of multiple filamentation. It is shown that the growth rate of this instability decreases when it is considered that the index of refraction has a stochastic behavior. A good qualitative agreement between the analytical and numerical results is obtained.

## 1. Introduction

A better understanding of the propagation of very powerful short laser pulses is a crucial issue for its wide range of applications in the remote sensing of chemical and biological agents as well as in directed energy applications, induced electrical discharges, lightning protection, long-range propagation of light bullets [1–9] ... For instance, the LIDAR (light detection and ranging) technique is an effective tool to investigate the atmosphere pollutants.

At high laser power, where laser power exceeds the critical power $P \gg P_{cr}$, the input beam, due to unavoidable beam irregularities, breaks up into $N$ filaments, with $N \sim P/P_{cr}$, each of them carrying about the critical power [10,11], Chaotic bundles of filaments are formed, and an erratic backscattering signal is observed. This multiple filamentation can play a positive or negative role depending on what we want to measure or the application we are targeting. Lightning discharge control requires long homogeneous plasma channels which are produced by the filament [12–15]. The detected backscattered nitrogen fluorescence from inside the filaments yields irregular changes from shot to shot which cannot be explained by fluctuation resulting from the initial laser pulse itself [16]. This irregularity is not only because of the difference in signal profiles [16–18]. Multiple filamentation was also studied to explain the fading of

the backscattering signal during periodically pulsed probing the atmosphere [17]. Thus, for many applications, filamentation control is an important issue and there exists a challenge in controlling the modulational instability which gives rise to multiple filamentation.

Focusing on low laser intensities, we show that the growth rate of the Bespalov and Talanov instability [19] can be reduced by turbulence. In this article, turbulence is outlined by introducing a stochastic refractive index in the nonlinear Schrödinger equation. A stochastic differential equation is solved considering refractive index fluctuations. A rough approximation used to average the Green's function of the set of equations is legitimized by deriving, in one case, the same result without using this approximation (Appendix A).

In this way, we show that the growth rate of the Bespalov and Talanov instability can be minimized by turbulence. This will prevent the buildup of multiple filamentation. Our assumption is that at the start of the beam, there are certain maxima of energy density (Fig. 1). These irregularities can be amplified by the modulational instability. Due to Kerr law, that is to say to the nonlinear part $\bar{n}_2I$ of the index of refraction [1], the light rays are confined along the positive refractive index gradients and some filaments are created [8,20,21].

The purpose of this article is to emphasize that multiple filamentation can be minimized when the atmosphere is turbulent through

<footnote><sup>a</sup> Department of Physics and Astronomy, The University of New Mexico, Albuquerque, NM, USA</footnote>

<footnote><sup>b</sup> Laboratoire de Spectroscopie Atomique, Moléculaire et Applications (LSAMA), Faculty of Sciences of Tunis, University of Tunis El Manar, Tunis, Tunisia</footnote>

<footnote><sup>c</sup> Laboratoire de l'atmosphère et des cyclones (LACy), CNRS, UMR 8105, Université de la Réunion, France</footnote>

<footnote>* Corresponding author.
E-mail address: alain.bourdier@gmail.com (A. Bourdier).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Multiple filamentation is initiated by the modulational instability.</figcaption></figure>

an analytical approach. This leads to a beam of light being more uniform, meaning it is almost free from intensity discontinuities.

## 2. Model equations

We consider an electromagnetic field propagating in the atmosphere. Its electric field has the following form

$$E[r=(x,y,z),t] = \psi(r,t) \expi(k_0 z - \omega_0 t) \frac{\hat{\epsilon}_z}{2} + cc, \quad (1)$$

where $\psi(r,t)$ varies slowly in time and $z$ and with $k_0 = n_0(\omega_0/c)$, where $n_0$ is the space average of the linear refractive index and $r = (x,y,z)$ the spatial dependence. We have $n = n_0 + 2n_2\langle E^2(r,t) \rangle + \delta n(r)$ where the angular brackets represent a time average and $\delta n(r)$ is a random function which characterizes the atmosphere turbulence, the nonlinear index which quantifies Kerr effect is considered here, $n_2$ is the second-order index of refraction [1,22–24]. The refractive index is assumed to depend on space only as a ray goes through the atmosphere so quickly that it does not change significantly during the crossing.

Fourier transforming the wave equation gives

$$\frac{\partial^2}{\partial z^2} E(r, \omega) + \Delta_{\perp} E(r, \omega) + \frac{\varepsilon(\omega)}{\varepsilon_0} \frac{\omega^2}{c^2} E(r, \omega) = 0, \quad (2)$$

the quantity $\varepsilon$ is the atmosphere dielectric constant and $\varepsilon_0$ the permittivity of free space. As the amplitude of the wave is assumed to vary slowly in time, $\psi$ contains no high-frequency components and we consider that $E(r, \omega) \approx (1/2)\psi(r, \omega - \omega_0)e^{ik_0 z}$. Moreover, the slowly varying amplitude in space approximation is made and $\partial^2\psi/\partial z^2$ can be dropped [22]. The wave equation, in terms of $\psi$, is

$$2ik_0 \frac{\partial}{\partial z} \psi + \Delta_{\perp} \psi + (k^2 - k_0^2) \psi = 0, \quad (3)$$

where $\psi = \psi(r, \omega - \omega_0)$ and $k(\omega) = \sqrt{\varepsilon(\omega)/\varepsilon_0\omega/c}$. As $\langle E^2(r, t) \rangle = (1/2)\psi\psi^* = (1/2)|\psi|^2$, we have now $n = n_0 + n_2|\psi|^2 + \delta n(r)$.

Then, the equation of propagation is reformulated in the time domain by using the inverse Fourier transform. Next, a retarded time is introduced. Assuming that the group velocity dispersion can be neglected, the equation for the propagation of the pulse in a turbulent atmosphere reads [22,25–27]

$$2ik_0 \psi_z[r=(x,y,z)] + \Delta_{\perp} \psi(r) + 2k_0^2 \frac{n_2}{n_0} |\psi|^2 \psi(r) = -2 \frac{k_0^2}{n_0} \delta n \psi(r). \quad (4)$$

To predict the amplification of irregularities initially present on the laser wavefront, a radial perturbation of $\psi$ is introduced [19,22]

$$\psi(r) = \psi_0(z) + a_1(z)e^{ik_{\perp} \cdot r} + a_{-1}(z)e^{-ik_{\perp} \cdot r}, \quad (5)$$

where $a_1(z)$ and $a_{-1}(z)$ are first-order quantities and $\pm k_{\perp}$ are the transverse components of the wavevector of the off-axis modes. We begin by assuming that fluctuations in the refractive index are in the form

$$\delta n(z) = \delta n_0 e^{i\phi(z)}, \quad (6)$$

where $\delta n_0$ is a first-order quantity and $\phi$ a stochastic quantity.

To find a solution to the wave equation [Eq. (4)], we first assume that $\psi_0$ satisfies

$$\frac{d\psi_0(z)}{dz} - i \frac{\omega_0}{c} \psi_0 (n_2 |\psi_0|^2 + \delta n_0 e^{i\phi}) = 0, \quad (7)$$

The solution is $\psi_0(z) = \bar{\psi} \exp \left[ i \left( \gamma z + \frac{\omega_0}{c} \delta n_0 \int_0^z e^{i\phi} dz \right) \right]$ where $\gamma = \frac{\omega_0}{c} n_2 \bar{\psi}^2$ and $\bar{\psi}$ is a real constant which is a zero-order quantity. We consider $\xi = \delta n_0 / (n_2 \bar{\psi}^2)$ is small and neglect the stochastic integral.

Always with the aim of satisfying the wave equation, the expansion for $\psi(r)$ is considered. The propagation equation is satisfied by setting to zero the terms in $e^{ik_{\perp} \cdot r}$ and in $e^{-ik_{\perp} \cdot r}$

$$\begin{aligned} 2ik_0 \frac{\partial a_1}{\partial z} - k_{\perp}^2 a_1 &= -k_0^2 \frac{2n_2}{n_0} \bar{\psi}^2 (2a_1 + a_{-1}^* e^{2i\gamma z}) - \frac{2k_0^2}{n_0} \delta n_0 e^{i\phi} e^{i\gamma z} a_1, \\ 2ik_0 \frac{\partial a_{-1}}{\partial z} - k_{\perp}^2 a_{-1} &= -k_0^2 \frac{2n_2}{n_0} \bar{\psi}^2 (2a_{-1} + a_1^* e^{2i\gamma z}) - \frac{2k_0^2}{n_0} \delta n_0 e^{i\phi} e^{i\gamma z} a_{-1}. \end{aligned} \quad (8)$$

Letting $a_{\pm 1} = \tilde{a}_{\pm 1} e^{i\gamma z}$, $\beta = k_{\perp}^2/2k_0$, $\tilde{\delta}n_0 = \delta n_0 e^{i\gamma z}$, and $\chi = k_0 \tilde{\delta}n_0/n_0$, the following form can be used to write this set of equations

$$\frac{d}{dz} \left( \frac{\tilde{a}_1}{\tilde{a}_{-1}^*} \right) = M_{\psi} \left( \frac{\tilde{a}_1}{\tilde{a}_{-1}^*} \right), \quad (9)$$

with

$$M_{\psi} = \begin{bmatrix} i(\gamma - \beta + \chi e^{i\phi}) & i\gamma \\ -i\gamma & -i(\gamma - \beta + \chi e^{-i\phi}) \end{bmatrix}. \quad (10)$$

At low intensity $\gamma$ is small and $\tilde{\delta}n_0$ does not vary significantly with $z$. In this case, $\chi$ may be supposed to be a constant. The stochastic phase $\phi$ is assumed to undergo the Kubo-Anderson process [28,29] which is set as a stepwise constant random function that jumps at randomly chosen distances between random step-values (Fig. 2).

The jumping positions of the stepwise constant random function $\phi(z)$ are uniformly and independently Poisson distributed with mean rate $\Delta\kappa$. We have $\langle e^{i\phi(z)} \rangle = 0$ and $e^{[i\phi(z)-i\phi(z)]} = e^{(-\Delta\kappa|z-z|)}$

## 3. Calculation of the growth rate

It is convenient to introduce the Green's function $G$ satisfying the matrix equation [28,29]

<figure><img src="image_5.png" /><figcaption>Fig. 2. Evolution of the phase $\phi$ versus the propagation distance.</figcaption></figure>
$$ \frac{dG(z, z')}{dz} = M_{\varphi}(z)G(z, z'), \qquad (11) $$

with $G(z', z) = I$, where $I$ is the identity matrix. Then, in terms of $G$, the solution of Eq. (9) reads

$$ A(z) = G(z, 0)A(0), \qquad (12) $$

with

$$ A(z) = \begin{pmatrix} \tilde{a}_1 \\ \tilde{a}_{-1}^* \end{pmatrix}. \qquad (13) $$

Based upon the Green's function equation we have $dG(z, 0)/dz = M_{\varphi}(z)G(z, 0)$. Then, we calculate the mean Green's function $\langle G(z, 0) \rangle$ on the different jumping positions $z_i$ and phases $\varphi_i$. Remembering that $\langle A(z) \rangle = \langle G(z, 0) \rangle A(0)$, we look for a master equation for this function. It can be written as the sum of two components: the first representing the part containing no jump between (0, z) and the second one integrating jumps occurring between (0, z).

Considering a Poisson distribution, the probability of having no jump between 0 and z is $\exp(-\Delta\kappa z)$. Then, the contribution to $\langle G(z, 0) \rangle$ may be written: $\langle G(z, 0) \rangle_{\text{nojump}} = \exp(-\Delta\kappa z)\langle G_0(z, 0) \rangle$, where $G_0(z, 0) = \exp[M_{\varphi}z]$ is obtained by solving Eq. (11) when $\varphi$ is assumed to be a constant and $\langle G_0(z, 0) \rangle = \langle G_0(z, 0) \rangle_{\varphi}$ is obtained by averaging over $\varphi$.

In the case when we have jumps, the jumping positions between 0 and z take place at the points $z_1 < z_2 < \dots < z_n$ and the semigroup property of the Green's function leads to

$$ G(z, 0) = G_0(z, z_n)G_0(z_n, z_{n-1})\dots G_0(z_1, 0), \qquad (14) $$

for a single realization. If $z'$ is the last jump position before z, it can be written

$$ G(z, 0) = G_0(z, z')G(z', 0), \qquad (15) $$

the probability for this last jump to occur between $z'$ and $z' + \delta z'$ is $\Delta\kappa\exp[-\Delta\kappa(z - z')]dz'$. Thus, the contribution of all the jumps to $\langle G(z, 0) \rangle$ is given by

$$ \begin{aligned} \langle G(z, 0) \rangle_{\text{jump}} &= \int_0^z \Delta\kappa dz' \exp(-\Delta\kappa(z - z')) \langle G_0(z, z') G(z', 0) \rangle_z \\ &= \int_0^z \Delta\kappa dz' \exp(-\Delta\kappa(z - z')) \langle G_0(z, z') \rangle \langle G(z', 0) \rangle_z, \end{aligned} \qquad (16) $$

as $G_0(z, z')$ and $G(z', 0)$ are independent. The symbol $\langle \cdot \rangle_z$ in Eq. (16) represents the conditional average of $A$ as a last jump occurred at $z'$. Moreover, as we consider a stationary problem, we have $G(z, z') = G(z - z', 0) = G(z - z')$. Thus, the average total Green's function is given by

$$ \langle G(z, 0) \rangle = \exp(-\Delta\kappa z) \langle G_0(z, 0) \rangle + \int_0^z \Delta\kappa dz' \exp( -\Delta\kappa(z - z')) \langle G_0(z - z') \rangle \langle G(z', 0) \rangle_z. \qquad (17) $$

We have $\langle G(z', 0) \rangle_z = \langle G(z', 0) \rangle$ as the process is independent on future jumps.

This equation is solved by performing a Laplace transform: $G(p) = \int_0^z dz e^{-p z} \langle G(z, 0) \rangle$. When considering the first term in the right-hand side of (17) which corresponds to the no jump case, we have

$$ G_0(p) = \int_0^z dz e^{-p z} \langle G(z, 0) \rangle_{\text{nojump}} = \int_0^z dz e^{-p z} \frac{1}{2\pi} \int_0^{2\pi} d\varphi \exp(M_{\varphi} z). \qquad (18) $$

Note that when we take the Laplace transform of a matrix equation of the form $\dot{G}_0(z, 0) = M_{\varphi}G_0(z, 0)$, we obtain $p\tilde{G}_0(p, 0) - G_0(0, 0) = M\tilde{G}_0(p, 0)$, thus $\tilde{G}_0(p, 0) = (pI - M_{\varphi})^{-1}$. Therefore, Eq. (18) reads

$$ G_0(p) = \langle \tilde{G}_0(p, 0) \rangle = \langle (pI - M_{\varphi})^{-1} \rangle. \qquad (19) $$

The Laplace transform of Eq. (17) gives

$$ \begin{aligned} G(p) &= [I - \Delta\kappa G_0(p + \Delta\kappa)]^{-1} G_0(p + \Delta\kappa) \\ &= [I - \Delta\kappa \langle ((p + \Delta\kappa)I - M_{\varphi})^{-1} \rangle]^{-1} \langle ((p + \Delta\kappa)I - M_{\varphi})^{-1} \rangle \\ &= N_1 N_2, \end{aligned} \qquad (20) $$

with $N_1 = [I - \Delta\kappa \langle ((p + \Delta\kappa)I - M_{\varphi})^{-1} \rangle]^{-1}$ and $N_2 = \langle ((p + \Delta\kappa)I - M_{\varphi})^{-1} \rangle$.

Finally, the growth rate of the perturbation will be derived by taking the inverse Laplace transformation of $G(p): \langle G(z, 0) \rangle =$

$$ \frac{1}{2\pi i} \int_{\Gamma - i\infty}^{\Gamma + i\infty} G(p) e^{p z} dp $$

### 3.1. Turbulence is not considered

It is assumed here that $\delta n_0 = 0$, then

$$ M_{\varphi} = M_0 = \begin{bmatrix} i(\gamma - \beta) & i\gamma \\ -i\gamma & -i(\gamma - \beta) \end{bmatrix}. \qquad (21) $$

Considering $K = \begin{pmatrix} a & b \\ c & d \end{pmatrix}$, then $K^{-1} = \frac{1}{ad-bc} \begin{pmatrix} d & -b \\ -c & a \end{pmatrix}$, consequently

$$ G_0(p) = (pI - M_0)^{-1} = \frac{1}{D_0} \begin{bmatrix} p + i(\gamma - \beta) & i\gamma \\ -i\gamma & p - i(\gamma - \beta) \end{bmatrix}, \qquad (22) $$

with $D_0 = p^2 + (\gamma - \beta)^2 - \gamma^2$. The poles of $G(p)$ are the zeros of $D_0$

$$ p = \pm p_0 = \pm \sqrt{\beta(2\gamma - \beta)}. \qquad (23) $$

The growth rate of the perturbation is derived by taking the inverse Laplace transformation of $G_0(p)$, $\langle G_0(z, 0) \rangle = G_0(z, 0) = \int_{\Gamma - i\infty}^{\Gamma + i\infty} \frac{1}{2\pi i} \int_{\Gamma - i\infty}^{\Gamma + i\infty} G_0(p) e^{p z} dp \, dz$. Applying the residue theorem: $G_0(z, 0) = \text{Res}(p_0) + \text{Res}(-p_0)$ where $\text{Res}(p_0)$ is the residue of $G_0(p)$ at $p_0$, only one pole gives a growth of the perturbation with $z$, for $z \gg 1$, we obtain

$$ G_0(z, 0) \sim \text{exp}(p_0 z), \qquad (24a) $$

thus

$$ A(z) = (\text{exp} p_0 z) A(0). \qquad (24b) $$

We have found again the result previously derived by Bespalov and Talanov [19]. Appendix B provides other ways to obtain this result.

### 3.2. Conclusions for the non-turbulent situation

The positive gain of the modulational instability is $p_0$ as far as $2\gamma - \beta > 0$, the value already derived by Bespalov and Talanov [19]. When $2\gamma - \beta < 0$, $p_0$ becomes imaginary which means that the energy transfer from the pump to the signals $a_1(z)e^{i\pm r}$ and $a_{-1}(z)e^{-i\pm r}$ is completely inhibited.

Figs. (4–8) show the good agreement between the numerical integration of Eq. (9) and the analytical result; the blue dashed line and the blue solid line one have the same slope.

This modulational instability can be studied in the framework of the four-wave mixing [22,30,31]. Introducing $\Delta k = 2(\gamma - \beta)$ the phase mismatch between the two signals and the cross-coupled polarization [22,31,32], the growth rate can be put in the following form $p_0 = \sqrt{\gamma^2 - \Delta k^2/4}$ [Appendix C]. When the mismatch $\Delta k$ is large enough the
growth rate becomes imaginary, and the instability is suppressed. When $\Delta\kappa = 0$, the growth rate is maximum, the four-wave mixing-process is phase matched [31]. Considering a stochastic index implies that the different wave vectors are stochastic, it can be considered that this situation is the superposition of several realizations which implies that the matching condition is not satisfied in many cases. As a result, the growth rate should be reduced.

### 3.3. Turbulence is considered

Let us verify that the growth rate of the instability is really affected by the stochastic index. We have

$$ (p'I - M_{\varphi})^{-1} = \frac{1}{D} \left[ \begin{matrix} p' + i(\gamma - \beta + \chi e^{-i\varphi}) & i\gamma \\ -i\gamma & p' - i(\gamma - \beta + \chi e^{i\varphi}) \end{matrix} \right], \quad (25) $$

with $p' = p + \Delta\kappa$ and $D = p'^2 + ip'\chi(e^{-i\varphi} - e^{i\varphi}) + (\gamma - \beta)^2 + (\gamma - \beta)\chi(e^{-i\varphi} + e^{i\varphi}) + \chi^2 - \gamma^2$.

Then, we have used the following rough approximation, $\langle x_i / y_i \rangle = [N(x_i) + \sum_{i=1}^N \delta x_i] / [N(y_i) + \sum_{i=1}^N \delta y_i] \approx \langle x_i \rangle / \langle y_i \rangle$ where N is the number of terms to average, $\delta x_i$ and $\delta y_i$ are fluctuations around the average. It is shown in Appendix A that when considering another form for refractive index fluctuation $\delta n$, the results obtained by using this approximation and the one given by an exact calculation are the same showing that this approximation is valid. It leads to

$$ N_2 = \langle (p'I - M_{\varphi})^{-1} \rangle = \frac{1}{\langle D \rangle} = \frac{1}{(p'^2 + (\gamma - \beta)^2 + \chi^2 - \gamma^2)} \left[ \begin{matrix} p' + i(\gamma - \beta) & i\gamma \\ -i\gamma & p' - i(\gamma - \beta) \end{matrix} \right], \quad (26) $$

and

$$ N_1 = [I - \Delta\kappa \langle (p'I - M_{\varphi})^{-1} \rangle]^{-1} = \frac{1}{R} \left[ \begin{matrix} 1 - [p' - i(\gamma - \beta)] \frac{\Delta\kappa}{\langle D \rangle} & i\gamma \frac{\Delta\kappa}{\langle D \rangle} \\ -i\gamma \frac{\Delta\kappa}{\langle D \rangle} & 1 - [p' + i(\gamma - \beta)] \frac{\Delta\kappa}{\langle D \rangle} \end{matrix} \right], \quad (27) $$

with $R = [1 - [p' - i(\gamma - \beta)] \frac{\Delta\kappa}{\langle D \rangle}] [1 - [p' + i(\gamma - \beta)] \frac{\Delta\kappa}{\langle D \rangle}] - (\gamma \frac{\Delta\kappa}{\langle D \rangle})^2$.

Thus, considering a small value of $\Delta\kappa$ and a large enough value of $\langle D \rangle$ so that $\Delta\kappa/\langle D \rangle$ can be assumed to be a small quantity, then, when $p \approx p_0$, we have $R \approx 1 - 2p'(\Delta\kappa/\langle D \rangle)$. If $\langle D \rangle$ is very small, then $R \propto 1/\langle D \rangle$. The simple poles of $G(p) = N_1 N_2$ [Eq. (20)] are the zeros of $\langle D \rangle$. If p is frankly different from $p_0$, the poles no longer exist and the growth rate fades. Finally, the growth rate modified by turbulence reads

$$ p = p_0 \left( 1 - \frac{1}{2} \frac{\chi^2}{p_0^2} \right) - \Delta\kappa. \quad (28) $$

The growth rate of the modulational instability is lessened by turbulence. It is significantly attenuated when $\chi \approx p_0$ and when $\Delta\kappa \approx p_0$.

Eqs. (9) were also integrated numerically in some physical circumstances considering that the evolution of the phase $\varphi$ versus the propagation is close to the one of the theoretical part (Fig. 1).

It is assumed that between two jumps phase value $\varphi$ is constant. The phase variations are described in Fig. 3. The values of $\varphi$ is determined by a random number and is supposed in a first step, to jump every constant distance $azk = 1/\Delta\kappa$.

For the sake of simplicity, the phase versus z was first assumed to be a staircase function y multiplied by $2\pi$, $\varphi$ is defined as follows: $y = -1$ when $-azk < z < 0$ and $y = 1$ when $0 < z < azk$.

The Fourier series for this function is

<figure><img src="image_6.png" /><figcaption>Fig. 3. Evolution of the phase $\varphi$ versus the propagation distance.</figcaption></figure>

$$ \varphi(z) = 8 \left[ \sin(\Delta\kappa\pi z) + \frac{1}{3}\sin(3\Delta\kappa\pi z) + \frac{1}{5}\sin(5\Delta\kappa\pi z) + \dots \right] \frac{1}{2n+1} \sin[(2n+1)\Delta\kappa\pi z] + \dots. \quad (29) $$

A first step was to approach the actual phase variations by considering many terms of Fourier's expansion.

Then the phase $\varphi$ was assumed to be a stochastic function. A stochastic fix to the Runge Kutta was implemented. Here, the phase $\varphi$ was given by $r \times 2\pi$ where r is a random number which should be drawn from the distribution $N(0, h)$ with zero mean and a variance of h the integration step. The deterministic dynamics was updated every integration step while the stochastic terms are updated using the square root of h [33–36]. We have $N(0, h) = \sqrt{h}N(0, 1)$, the fix was therefore obtained by multiplying the stochastic phase by the square root of the space step.

To validate the good agreement between the analytical and the numerical results different parameters were considered. First, the case of a CO2 laser propagating in the atmosphere was considered. In these simulations, we took, $\lambda = 10.6\mu m$, $\beta = \gamma$ and $n_2 = 1.43 \times 10^{-26}m^2/V^2$ [37]. Then, the case of a long pulse neodymium laser ($\lambda = 1.06\mu m$) was considered with $I = 2 \times 10^{10}W/cm^2$ [38] and $n_2 = 1.035 \times 10^{-25}m^2/V^2$ [8,39,40].

Figs. (4–8) show the evolution of the dimensionless quantity $\hat{I} = |\tilde{a}_1|$ as a function of z calculated in different ways. The gains p and $p_0$ are our two analytical results, $\exp(p_0z)$ is represented by a blue dashed line which shows the growth rates of the instability when no turbulence is considered, when turbulence is taken into account, a red dashed line is used to represent the quantity $\exp(pz)$. The normalized amplitude of the instability $|\tilde{a}_1|$ is also obtained by integrating numerically Eq. (9) with a simple fourth order Runge Kutta. The stochastic fix is used in the stochastic case, the results are represented by a solid red line. In the nonturbulent case they are represented by a solid blue line. Different runs with different phase $\varphi$-distributions were considered, the results obtained are always similar to those displayed in these figures. The black solid lines are obtained with different number of terms in the Fourier expansion [Eq. (29)].

Fig. (4a) shows the results obtained when $\chi$ is a function of z, while Fig. 4b shows the results obtained under the assumption that it is constant. We began to suppose that the distance azk between two jumps is constant.

It is apparent that considering $\chi$ is a constant makes sense, at least, over a distance close to the half-period of $\cos\gamma z$. In these cases, the growth obtained for the instability cannot be compared to the analytical result [Eq. (28)] as $\Delta\kappa$ is too large.

These last results show that in a certain domain, the growth rate of the modulational instability is indeed reduced by the turbulence of the atmosphere. We will now ensure that our numerical results are relevant by comparing them to an analytical result in the following.

From this point, $\chi$ is a constant, the values of the various parameters examined do not strictly match the frame we defined in our hypotheses.
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** $\lambda = 10.6\mu m$, $I = 1.5 \times 10^{11}W/cm^2$, $\delta n_0 / n_0 = 7 \times 10^{-10}$ and azk = 20, $\gamma = 9.5 \times 10^{-3}$, $\xi = 4.36 \times 10^{-2}$.</figcaption></figure>

The blue dashed line is the analytic result when no turbulence is considered $[\exp(p_0z)]$.

The black solid lines are obtained with different number of terms in the Fourier expansion.

The blue solid line is obtained by integrating numerically Eq. (9) when there is no turbulence ($|\tilde{a}_1|$). The red solid line is obtained when there is turbulence with a simple fourth order Runge Kutta by using a stochastic fix ($|\tilde{a}_1|_{\text{FIX}}^T$).

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** $\lambda = 10.6\mu m$, $I = 1.5 \times 10^{12}W/cm^2$, $\delta n_0 / n_0 = 3.3 \times 10^{-9}$, azk = 30.</figcaption></figure>

$\gamma = 9.58 \times 10^{-2}$, $\xi = 2 \times 10^{-2}$.

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** $\lambda = 10.6\mu m I = 1.5 \times 10^{12}W/cm^2$, $\delta n_0 / n_0 = 3 \times 10^{-9}$, azk = 40. $\gamma = 9.58 \times 10^{-2}$, $\xi = 1.85 \times 10^{-2}$.</figcaption></figure>

The brown curve is obtained with a simple fourth order Runge Kutta by using a stochastic fix ($|\tilde{a}_1|_{\text{FIX}}^T$) and a constant azk (azk=40). The red solid line is obtained in the same way but with a random azk ($|\tilde{a}_1|_{\text{FIX},1/\text{azkrandom}}^T$).

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** $\lambda = 1.06\mu m I = 2 \times 10^{10}W/cm^2$, $\delta n_0 / n_0 = 3.5 \times 10^{-10}$, azk = 40. $\gamma = 9.8 \times 10^{-2}$, $\xi = 2.24 \times 10^{-2}$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** $\lambda = 1.06\mu m I = 2 \times 10^{10}W/cm^2$, $\delta n_0 / n_0 = 3.5 \times 10^{-10}$, azk = 25. $\gamma = 1.13 \times 10^{-1}$, $\xi = 1.9 \times 10^{-2}$.</figcaption></figure>
Our goal is to ensure that the analytical and numerical solutions for our set of equations are consistent [Eq. (9)].

In Fig. 5, results corresponding to high values of z and high values of $|\tilde{a}_1|$ are shown.

A case when the distance between two jumps of the phase is constant is compared to a case when it is stochastic is shown in Fig. 6. The stochastic distance azk is obtained by multiplying the constant value by a random number between 0 and 1. The integrations of the stochastic equation in the two cases give similar results.

For $\lambda = 10.6\mu m$, in these cases, the different numerical integrations of Eq. (9) confirm that the growth of the instability is reduced by turbulence. Our overall results indicate that there is good qualitative agreement between the analytical approach and the numerical approach.

Let us now move on to a value of the wavelength more accessible to an experiment: $\lambda = 1.06\mu m$. Fig. 7 shows results when turbulence is considered and when azk is a constant

A case when the distance between two jumps of the phase is constant was also compared to one when it is stochastic. The results are similar to those shown in Fig. 6.

Thus, in the cases when $\lambda = 1.06\mu m$, our results still show a good qualitative agreement between the analytical approach and the numerical one which was obtained using the same fix as before. We find that the growth of the instability is still reduced by turbulence.

It was also assumed that the refractive index is in the form

$$ \delta n(z) = \delta n_0 \cos\varphi(z), \quad (30) $$

contrary to a previous situation, the refractive index has no imaginary part.

To find a solution to the wave equation [Eq. (4)], we first assume that $\psi_0$ satisfies

$$ \frac{d\psi_0(z)}{dz} - i \frac{\omega_0}{c} \psi_0 (n_2 |\psi_0|^2 + \delta n_0 \cos\varphi) = 0. \quad (31) $$

The solution is $\psi_0(z) = \bar{\psi} \exp \left[ i \left( \gamma z + \frac{\omega_0}{c} \delta n_0 \int_0^z \cos\varphi dz \right) \right]$. We have also considered $\xi = \delta n_0 / n_2 \bar{\psi}^2$ is small and neglected the stochastic integral.

The expansion for $\psi(r)$ is considered again assuming that $\chi$ varies slowly with $z$, the propagation equation is satisfied by setting to zero the terms in $e^{ik_{\perp} \cdot r}$ and in $e^{-ik_{\perp} \cdot r}$. Then, introducing the same parameters and variables as before, the wave equation for the off-axis modes is obtained and reads

$$ \frac{d}{dz} \begin{pmatrix} \tilde{a}_1 \\ \tilde{a}_{-1}^* \end{pmatrix} = M_{\varphi} \begin{pmatrix} \tilde{a}_1 \\ \tilde{a}_{-1}^* \end{pmatrix}, \quad (32) $$

with

$$ M_{\varphi} = \begin{bmatrix} i(\gamma - \beta + \chi \cos\varphi) & i\gamma \\ -i\gamma & -i(\gamma - \beta + \chi \cos\varphi) \end{bmatrix}. \quad (33) $$

As $\delta\tilde{n}_0$ does not vary significantly with $z$, $\chi$ still may be supposed to be a constant. The stochastic phase $\varphi$ is still assumed to undergo the Kubo-Anderson process.

The following growth can be found by carrying out the same analytical work as before

$$ p = p_0 \left( 1 - \frac{1}{4} \frac{\chi^2}{p_0^2} \right) - \Delta\kappa, \quad (34) $$

thus, the growth rate of modulational instability is indeed reduced by turbulence. Now moving to the numerical approach, considering $\lambda = 1.06\mu m$, $I = 2 \times 10^{10} \text{W/cm}^2$ and $n_2 = 1.2 \times 10^{-25} \text{m}^2/\text{V}^2$ a numerical integration of Eq. (32) was achieved. Fig. 8 shows results when turbulence is considered and when the distance between two jumps, azk, is a constant. It also shows that with real refractive index fluctuations, the effect of the turbulence is still to decrease the growth of the modulational instability.

To summarize, the results presented above exhibit some consistency between our analytical and numerical results. For low intensities, they show that the modulational instability originally described by V.I. Bespalov and V.I. Talanov [19] is attenuated by turbulence, which is consistent with other studies [20].

## 4. Conclusions

The impact of turbulence on the propagation of a laser pulse in the atmosphere, multiple filamentation have already been studied extensively [9,20,25,41–43]. Here, we have studied its influence of turbulence on the multiple filamentation process when a laser beam moves through the atmosphere. Its action on the modulational instability of a plane wave was investigated within the framework of the nonlinear Schrödinger equation. First, we have established the wave propagation equation when considering a stochastic refractive index perturbation to simulate turbulence. Our focus was on low laser intensities. The stochastic index was expressed by using the Kubo-Anderson process. The Green’s function of the wave equation was introduced. An equation for the average Green’s function was determined using its semigroup property. Thus, the growth rate of the instability which begets the laser beam filamentation was calculated considering different forms for the stochastic part of the refractive index, one of them allowing a rigorous method. This was used to justify the rough approximation made to determine the average Green’s function and its poles.

The stochastic wave equation has been solved numerically in a very close way first by replacing the stochastic part of the equation by a function giving a similar evolution. Then, a fix to the Runge Kutta method to solve the stochastic set of equations which describe the instability is implemented. This fix was obtained by multiplying the stochastic phase by the square root of the space step.

To test our numerical results, we have used the analytical solution.

A global agreement between analytical and numerical results was reached. We have shown analytically and numerically a trend: turbulence decreases the growth rate of the Bespalov and Talanov instability. Therefore, if we extrapolate the conclusion derived from our model and apply it to high powers, turbulence should cancel out most of the multiple filamentation phenomenon.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

No data was used for the research described in the article.
Appendix A

The fluctuations of the refractive index are defined in a different way

The refractive index fluctuations are assumed to be in the form

$$
\delta n(r) = \delta n_0 \left[ e^{i(k_{\perp} \cdot r + \varphi(z))} + e^{-i(k_{\perp} \cdot r + \varphi(z))} \right]. \quad (A1)
$$

This form was chosen to find an analytical solution in two different manners, one using the rough approximation and the other not. The first goal was to legitimate the rough approximation. The second objective was to continue examining the numerical resolution of our stochastic differential equations.

The calculations are conducted using the nonlinear Schrödinger equation with this expression for $\delta n(r)$ exactly as in paragraph 2.

The same expansion for $\psi(r)$ as in paragraph 2 is considered. The following equation must be verified

$$
\psi_0(z) = \bar{\psi} \exp(i\gamma z). \tag{A2}
$$

We consider that $n_2 \bar{\psi}^2 a_{\pm 1} \sim \delta n_0 \bar{\psi} \gg \delta n_0 a_{\pm 1}$ which means that $\delta n_0$ is supposed to be very small. Then, the propagation equation is satisfied by setting to
zero the terms in $e^{ik_{\perp} \cdot r}$ and in $e^{-ik_{\perp} \cdot r}$. We obtain

$$
\begin{align}
& 2ik_0 \frac{\partial a_1}{\partial z} - k_{\perp}^2 a_1 = -k_0^2 \frac{2n_2}{n_0} \bar{\psi}^2 (2a_1 + a_{-1}^* e^{2i\gamma z}) - \frac{2k_0^2}{n_0} \bar{\psi} \delta n_0 e^{i\gamma z} e^{i\varphi} a_1, \nonumber \\
& 2ik_0 \frac{\partial a_{-1}}{\partial z} - k_{\perp}^2 a_{-1} = -k_0^2 \frac{2n_2}{n_0} \bar{\psi}^2 (2a_{-1} + a_1^* e^{2i\gamma z}) - \frac{2k_0^2}{n_0} \bar{\psi} \delta n_0 e^{i\gamma z} e^{-i\varphi} a_{-1}. \tag{A3}
\end{align}
$$

Introducing the same parameters and variables as before, the matrix $M_{\varphi}$ becomes

$$
M_{\varphi} = \begin{bmatrix} i(\gamma - \beta + \chi e^{i\varphi}) & i\gamma \\ -i\gamma & -i(\gamma - \beta + \chi e^{i\varphi}) \end{bmatrix}, \quad (A4)
$$

with here $\chi = k_0 \bar{\psi} \delta \tilde{n}_0 / n_0$

1) The rough approximation is used

We find

$$
(pI - M_{\varphi})^{-1} = \frac{1}{D} \begin{bmatrix} p' + i(\gamma - \beta + \chi e^{i\varphi}) & i\gamma \\ -i\gamma & p' - i(\gamma - \beta + \chi e^{i\varphi}) \end{bmatrix} \quad (A5)
$$

Let us calculate $G(p) = N_1 \cdot N_2$, thus

$$
N_2 = \langle (p'I - M_\varphi)^{-1} \rangle = \frac{1}{(p'^2 + (\gamma - \beta)^2 - \gamma^2)} \begin{bmatrix} p' + i(\gamma - \beta) & i\gamma \\ -i\gamma & p' - i(\gamma - \beta) \end{bmatrix}, \quad (A6)
$$

and

$$
N_1 = [I - \Delta\kappa \langle (pI - M_\varphi)^{-1} \rangle]^{-1} = \frac{1}{R} \begin{bmatrix} 1 + [p' - i(\gamma - \beta)] \frac{\Delta\kappa}{\langle D \rangle} & i\gamma \frac{\Delta\kappa}{\langle D \rangle} \\ -i\gamma \frac{\Delta\kappa}{\langle D \rangle} & 1 - [p' + i(\gamma - \beta)] \frac{\Delta\kappa}{\langle D \rangle} \end{bmatrix}, \quad (A7)
$$

with $R = [1 + |p'| - i(\gamma - \beta)] \frac{\Delta\kappa}{\langle D \rangle} [1 - |p'| + i(\gamma - \beta)] \frac{\Delta\kappa}{\langle D \rangle} + (\gamma \frac{\Delta\kappa}{\langle D \rangle})^2$ and $\langle D \rangle = p'^2 + (\gamma - \beta)^2 - \gamma^2$.

As previously, it can be shown that the zeros of $\langle D \rangle$ are poles of $G(p)$. Thus, the growth rate modified by turbulence is given by

$$
p \approx p_0 - \Delta\kappa. \tag{A8}
$$

We still find that the growth rate of the modulational instability is attenuated by turbulence.

2) The growth rate is calculated again without using the rough approximation.

To discuss the validity of the rough approximation we must calculate the pole of ⟨G(p)⟩ and consequently the growth rate. Here, we shall calculate it in another way thanks to this form for the index perturbation. We start with Eq. (20), we express the two terms which when multiplied gives G(p). First, we have
$$N_2 = \langle (p'I - M_\varphi)^{-1} \rangle = p^{-1} \langle (I + p'^{-1} M_\varphi + p'^{-2} M_\varphi^2 + \dots + p'^{-n} M_\varphi^n + \dots) \rangle, \quad (A9)$$

with

$$\langle M_\varphi \rangle = M_0 = \begin{bmatrix} i(\gamma - \beta) & i\gamma \\ -i\gamma & -i(\gamma - \beta) \end{bmatrix}, \quad (A10)$$

$$M_{\varphi}^{2} = [\gamma^{2} - (\gamma - \beta + \chi e^{i\varphi})^{2}] I, \quad (A11)$$

and

$$M_{\varphi}^{3} = [\gamma^{2} - (\gamma - \beta + \chi e^{i\varphi})^{2}] M_{\varphi}. \quad (A12)$$

Consequently

$$\langle M_{\varphi}^{2} \rangle = [\gamma^{2} - (\gamma - \beta)^{2}] I, \quad (A13)$$

and

$$\langle M_{\varphi}^{3} \rangle = [\gamma^{2} - (\gamma - \beta)^{2}] M_{0}, \quad (A14)$$

we also have $\langle M_{\varphi}^{4} \rangle = [\gamma^{2} - (\gamma - \beta)^{2}]^{2} I$ and $\langle M_{\varphi}^{5} \rangle = [\gamma^{2} - (\gamma - \beta)^{2}]^{2} M_{0}$. Finally, we find

$$\begin{aligned} N_2 &= \langle (p'I - M_\varphi)^{-1} \rangle = p'^{-1} [I + p'^{-1} M_0 + p'^{-2} (\gamma^2 - (\gamma - \beta)^2) I + p'^{-3} [\gamma^2 - (\gamma - \beta)^2] M_0 + \dots] \\ &= p'^{-1} I [1 + p'^{-2} (\gamma^2 - (\gamma - \beta)^2) + p'^{-4} (\gamma^2 - (\gamma - \beta)^2)^2 + p'^{-6} (\gamma^2 - (\gamma - \beta)^2)^3 + \dots] \\ &\quad + p'^{-2} M_0 [1 + p'^{-2} (\gamma^2 - (\gamma - \beta)^2) + p'^{-4} (\gamma^2 - (\gamma - \beta)^2)^2 + p'^{-7} (\gamma^2 - (\gamma - \beta)^2)^3 + \dots] \\ &= \left[ \frac{p'^{-1}}{1 - p'^{-2} [\gamma^2 - (\gamma - \beta)^2]} \right] (p'^{-1} I + p'^{-2} M_0). \end{aligned} \quad (A15)$$

The second term necessary for the expression of $G(p)$ is

$$\begin{aligned} N_1 &= [I - \Delta\kappa \langle (p'I - M_\varphi)^{-1} \rangle]^{-1} \\ &= I + \Delta\kappa \langle (p'I - M_\varphi)^{-1} \rangle + [\Delta\kappa \langle (p'I - M_\varphi)^{-1} \rangle]^2 + \dots + [\Delta\kappa \langle (p'I - M_\varphi)^{-1} \rangle]^n + \dots \end{aligned} \quad (A16)$$

If $\Delta\kappa$ is small enough only the term $I$ can be considered. As $G(p) = N_2 \cdot N_1$, $p'_0 = \sqrt{\beta(2\gamma - \beta)}$ which is a pole of $\langle (p'I - M_\varphi)^{-1} \rangle$ is also pole of $G(p)$. It gives the growth rate of the instability

$$p'_{0} = \sqrt{\beta(2\gamma - \beta)}, \quad (A17)$$

that is to say

$$p \approx p_0 - \Delta\kappa, \quad (A18)$$

which is the previous result [Eq. (A8)]. It means that the rough approximation that we made previously to calculate the poles of $\langle G(p) \rangle$ leads to a good result.

Once again, we have compared our analytical results to the numerical results in the case of a $\lambda = 10.6\mu m$ wavelength laser pulse and in the case of a $\lambda = 1.06\mu m$ laser pulse. Here, we don't pay attention to the fact the different parameters fall outside or inside the framework defined by our assumptions. Our goal is just to emphasize the agreement between the analytical and numerical solutions of our equations.

The following results were obtained for $\lambda = 10.6\mu m$, $\delta n_0 / n_0 = 1 \times 10^{-16}$, $n_2 = 1.43 \times 10^{-26} m^2/V^2$, $I = 5 \times 10^{12} W/cm^2$ and $\text{azk} = 10$ when constant (Fig. A1)

The red solid line is obtained with a fourth order Runge Kutta by using a stochastic fix ($|\tilde{a}_1|_{\text{Fix},1/\Delta\text{krandom}}^T$) and a random $\text{azk}$. The brown solid one is obtained using a constant $\text{azk}(|\tilde{a}_1|_{\text{Fix}}^T)$.

In the case when $\lambda = 10.6\mu m$, $\delta n_0 / n_0 = 1 \times 10^{-16}$, $n_2 = 1.035 \times 10^{-25} m^2/V^2$, $I = 1 \times 10^{11} W/cm^2$ and $\text{azk} = 10$ we obtained (Fig. A2)

Figs. (A1 and A2) show that, in these cases, there is a good qualitative agreement between the analytical result and the numerical integration using the same fix as before. We find that the growth of the instability is still reduced by turbulence.
<figure><img src="image_12.png" /><figcaption>Fig. A1. $\lambda = 10.6\mu m$, $I = 5 \times 10^{12}W/cm^2$, $\delta n_0 / n_0 = 1 \times 10^{-16}$, $\text{azk} = 10$, $n_2 = 1.43 \times 10^{-26}m^2/V^2$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. A2. $\lambda = 1.06\mu m\sqrt{2}$, $I = 1 \times 10^{11}W/cm^2$, $\delta n_0 / n_0 = 1 \times 10^{-16}$, $\text{azk} = 10$, $n_2 = 1.035 \times 10^{-25}m^2/V^2$. The red solid line is obtained by using a stochastic fix ($|\tilde{a}_1|_{\text{Fix},1/\Delta\text{krandom}}^T$) and a random $\text{azk}$. The brown solid one is obtained using a constant $\text{azk}(|\tilde{a}_1|_{\text{Fix}}^T)$.</figcaption></figure>

## Appendix B

*Other ways to calculate the growth rate of the instability in the non-turbulent case*

In the non-turbulent case, we can expand $G_0(p)$ in the following way

$$G_0(p) = (pI - M_0)^{-1} = p^{-1}(I + p^{-1}M_0 + p^{-2}M_0^2 + \dots + p^{-n}M_0^n + \dots), \qquad (B1)$$

with

$$M_0 = \begin{bmatrix} i(\gamma - \beta) & i\gamma \\ -i\gamma & -i(\gamma - \beta) \end{bmatrix}. \qquad (B2)$$

We have

$$M_0^2 = [\gamma^2 - (\gamma - \beta)^2]I, \qquad (B3)$$

and

$$M_0^3 = [\gamma^2 - (\gamma - \beta)^2]M_0, \qquad (B4)$$

thus
$$
\begin{align}
G_0(p) &= p^{-1} \left\{ \sum_{n=0}^{\infty} \left[ (p^{-2}[\gamma^2 - (\gamma - \beta)^2])^n I + p^{-1} \sum_{n=0}^{\infty} (p^{-2}[\gamma^2 - (\gamma - \beta)^2])^n M_0 \right] \right\} \tag{B5} \\
&= \left[ \frac{p^{-1}}{1 - p^{-2}[\gamma^2 - (\gamma - \beta)^2]} \right] (I + p^{-1} M_0). \nonumber
\end{align}
$$

We have three poles: $p = 0$ and $p = \pm p_0 = \pm \sqrt{\beta(2\gamma - \beta)}$. Taking the inverse Laplace transformation, we obtain for $z \gg 1$

$$
G_0(z, 0) \sim \exp \pm (p_0 z) \tag{B6}
$$

Thus, the previous result obtained in the nonturbulent case in paragraph 3.1 has been derived again. The spirit of this method is important as it is applied in Appendix A.

Let us give another way to find this result. It can be easily shown that

$$
|\tilde{a}_1|^2 - |\tilde{a}_{-1}^*|^2 = q, \tag{B7}
$$

where $q$ is a constant. Then, we let $\tilde{a}_1 = q^{1/2}e^{i\varphi_1}\sinh\theta$ and $\tilde{a}_{-1}^* = q^{1/2}e^{i\varphi_2}\cosh\theta$. Here, $\varphi_1$ and $\varphi_2$ are two real constants. Assuming $\theta \gg 1$, considering the real part of Eq. (9) for $\tilde{a}_1$, we find

$$
\frac{d\theta}{dz} = -\gamma \sin(\varphi_2 - \varphi_1), \tag{B8}
$$

and the imaginary part leads to

$$
\gamma - \beta = -\gamma \cos(\varphi_2 - \varphi_1). \tag{B9}
$$

Combining (B8) and (B9) gives

$$
\theta = \pm p_0 z, \tag{B10}
$$

which means that the same result for the growth rate of the instability is found.

**Appendix C**

Influence of the phase mismatch on the growth rate of the modulational instability

The form of the laser electric field is assumed to be as follows

$$
E(r, t) = E(r)\exp(-i\omega_0 t) \tag{C1}
$$

The spatial dependence is the sum of three terms

$$
\begin{align}
E(r) &= E_0(r) + E_1(r) + E_{-1}(r) \tag{C2} \\
&= [A_0(z) + a_1(z)\exp(ik_\perp r) + a_{-1}(z)\exp(-ik_\perp r)]\exp(ik_0 z). \nonumber
\end{align}
$$

The wave equation reads

$$
\left( \Delta - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} \right) E(r, t) = \mu_0 \frac{\partial^2}{\partial t^2} (P_L + P_{NL}), \quad (C3)
$$

where $P_{NL}$ is the nonlinear polarization. We have: $P_{NL} = \epsilon_0\chi^{(3)}[E]^2 E$ where $\chi^{(3)}$ is the nonlinear susceptibility. One part of the polarization is phase matched to the high intensity part of the field, it is given by [22,30]

$$
P_{0NL} = \epsilon_0 \chi^{(3)} |E_0|^2 E_0, \tag{C4}
$$

and the part that is phase matched to the sidemodes is

$$
P_{\pm NL} = \epsilon_0 \chi^{(3)} [2|E_0|^2 E_{\pm 1} + E_0^2 E_{\pm 1}^*]. \tag{C5}
$$

We consider the pump wave, and the two signals are in the following form

$$
\begin{align}
E_0 &= A_0 \exp(ik_0 z - i\omega t), \notag \\
E_{\pm} &= a_{\pm} \exp(ik_{\pm} z - i\omega t). \tag{C6}
\end{align}
$$

When turbulence is not considered, Eqs (9) have a solution in the form.

$$
\begin{align}
\tilde{a}_1 &= \tilde{a}_1 \exp[i(\gamma - \beta)z], \notag \\
\tilde{a}_{-1}^* &= \tilde{a}_{-1} \exp[i(\gamma - \beta)z], \tag{C7}
\end{align}
$$

which must satisfy the following equations
$$ \frac{d^2}{dz^2} \hat{a}_{-1} + i\Delta k \frac{d}{dz} \hat{a}_{-1} - \gamma^2 \hat{a}_{-1} = 0, \qquad (C8) \\ \frac{d}{dz} \hat{a}_1 = i\gamma \hat{a}_{-1}. $$

The solution for the original variable is

$$ a_{\pm 1} \sim \exp \left[ \sqrt{\left( \gamma^2 - \frac{(\Delta k)^2}{4} \right) z} \right] \exp \left[ i \left( \frac{\Delta k}{2} + 2\gamma - \beta \right) z \right], \qquad (C9) $$

thus, we have: $k_{0T} = k_0 + \Delta k_{NL} = k_0 + \gamma$ and $k_{\pm} = k_0 + \frac{\Delta k}{2} - \beta + 2\gamma$.

The cross-coupled polarizations $E_0^2 E_{-1}^*$ and $E_0^2 E_{+1}^*$ must be phase matched to $E_{\mp}$. Thus, the mismatch must be

$$ \Delta \bar{k} = 2k_{0T} - 2k_{\pm} = -2\Delta k. \qquad (C10) $$

We find again that the growth rate is maximum is maximum when $\Delta k = 0$.

## References

[1] R.Y. Chiao, E. Garmire, C.H. Townes, Self-trapping of optical beams, Phys. Rev. Lett. 13 (1964) 479-482.

[2] A. Braun, G. Korn, X. Liu, D. Du, J. Squier, G. Mourou, Self-channeling of high-peak-power femtosecond laser pulses in air, Opt. Lett. 20 (1995) 73-75.

[3] P. Sprangle, E. Esarey, J. Krall, Self-guiding and stability of intense optical beams in gases undergoing ionization, Phys. Rev. E 54 (1996) 4211-4232.

[4] A. Brodeur, C.Y. Chien, F.A. Ilkov, S.L. Chin, O.G. Kosareva, V.P. Kandidov, Moving focus in the propagation of ultrashort laser pulses in air, Opt. Lett. 22 (1997) 304-306.

[5] S.L. Chin, S. Petit, F. Borne, K. Miyazaki, The White Light Supercontinuum Is Indeed an Ultrafast White Light Laser, JPN J. Appl. Phys. 38 (1999) L126-L128.

[6] J.R. Peñano, P. Sprangle, B. Hafizi, A. Ting, D.F. Gordon, C.A. Kapetanakos, Propagation of ultra-short, intense laser pulses in air, Phys. Plasmas 11 (2004) 2865-2874.

[7] A. Couairon, Light bullets from femtosecond filamentation, Eur. Phys. J. D 27 (2003) 159-167.

[8] A. Couairon, A. Mysyrowicz, Femtosecond filamentation in transparent media, Phys. Rep. 441 (2007) 47-189.

[9] J. Guo, L. Sun, J. Liu, B. Shang, S. Tao, N. Zhang, L. Lin, Z. Zhang, Beam Wander Restricted by Nonlinearity of Femtosecond Laser Filament in Air, Sensors 22 (2022) 4995.

[10] G. Méchain, A. Couairon, Y.-B. André, C. D'amico, M. Franco, B. Prade, S. Tzortzakis, A. Mysyrowicz, R. Sauerbrey, Long-range self-channeling of infrared laser pulses in air; a new propagation regime without ionization, Appl. Phys. B 79 (2004) 379-382.

[11] M. Rodriguez, R. Bourayou, G. Méjean, J. Kasparian, J. Yu, E. Salmon, A. Scholz, B. Stecklum, J. Eislöffel, U. Laux, A.P. Hatzes, R. Sauerbrey, L. Wöste, J.P. Wolf, Kilometer-range nonlinear propagation of femtosecond laser pulses, Phys. Rev. E 69 (2004), 036607.

[12] R. Ackermann, K. Stelmaszczyk, P. Rohwetter, G. Méjean, E. Salmon, J. Yu, J. Kaspariana, G. Méchain, V. Bergmann, S. Schaper, B. Weise, T. Kumm, K. Rethmeier, W. Kalkner, L. Wöste, J.P. Wolf, Triggering and guiding of megavolt discharges by laser-induced filaments under rain conditions, Appl. Phys. Lett. 85 (2004) 5781-5783.

[13] X.M. Zhao, J.-C. Diels, C.Y. Wang, J.M. Elizondo, Femtosecond Ultraviolet Laser Pulse Induced Lightning Discharges in Gases, IEEE J. Quant. Electron. 31 (1995) 599-612.

[14] F. Vidal, D. Comtois, C.-Y. Chien, A. Desparois, B. La Fontaine, T.W. Johnston, J.-C. Kieffer, H.P. Mercure, H. Pépin, F.A. Rizk, Modeling the Triggering of Streamers in Air by Ultrashort Laser Pulses, I.E.E.E. Trans, Plasma Sci. 28 (2000) 418-433.

[15] O.G. Kosareva, N.A. Panov, N. Akozbek, V.P. Kandidov, Q. Luo, S.A. Hosseini, W. Liu, J.-F. Gravel, G. Roy, S.L. Chin, Controlling a bunch of multiple filaments by means of a beam diameter, Appl. Phys. B 82 (2006) 111-122.

[16] S.A. Hosseini, Q. Luo, B. Ferland, W. Liu, S.L. Chin, O.G. Kosareva, N.A. Panov, N. Aközbek, V.P. Kandidov, Competition of multiple filaments during the propagation of intense femtosecond laser pulses, Phys Rev A 70 (2004), 033802.

[17] V.P. Kandidov, O.G. Kosareva, S.A. Shlenov, N.A. Panov, V.Yu. Fedorov, A. E. Dormidonov, Dynamic small-scale self-focusing of a femtosecond laser pulse, Quantum Elec. 35 (2005) 59-64.

[18] V.P. Kandidov, N. Akozbek, M. Scalora, O.G. Kosareva, A.V. Nyakk, Q. Luo, S. A. Hosseini, S.L. Chin, Towards a control of multiple filamentation by spatial regularization of a high-power femtosecond laser pulse, Appl. Phys. B 80 (2005) 267-275.

[19] V.I. Bespalov, V.I. Talanov, Filamentary structure of light beams in nonlinear liquids, JETP Lett 3 (1966) 307-310.

[20] V.P. Kandidov, S.A. Shlenov, O.G. Kosareva, Filamentation of high-power femtosecond laser radiation, Quantum Elec. 39 (2009) 205-228.

[21] Houard, A. 2015. Filamentation en regime TW et applications. Habilitation thesis, University Paris 6.

[22] R.W. Boyd, Nonlinear Optics, 2nd edition, Academic Press, Amsterdam, 2002.

[23] J.C. Diels, W. Rudolph, Ultrashort Laser Pulse Phenomena, Elsevier, 2006.

[24] N. Bloembergen, Nonlinear Optics, Benjamin W.A., Inc., New York, 1965.

[25] P.M. Lushnikov, N. Vladimirova, Toward defeating diffraction and randomness for laser beam propagation in yurbulent atmosphere, JETP Lett. 108 (2018) 571-576.

[26] J.M. Martin, S.M. Flatté, Intensity images and statistics from numerical simulation of wave propagation in 3-D random media, Appl. Opt. 27 (1988) 2111-2126.

[27] S.L. Chin, A. Talebpour, J. Yang, S. Petit, V.P. Kandidov, O.G. Kosareva, M. P. Tamarov, Filamentation of femtosecond laser pulses in turbulent air, Appl. Phys. B-lasers and optics 74 (2002) 67-76.

[28] A. Brissaud, U. Frisch, Solving linear stochastic differential equations, J. Math. Phys. 15 (1974) 524-534.

[29] J.J. Thomson, Finite-bandwidth effects on the parametric instability in an inhomogeneous plasma, Nucl. Fusion 15 (1975) 237-247.

[30] R.Y. Chiao, P.L. Kelley, E. Garmire, Stimulated four-photon interaction and its influence on stimulated Rayleigh-wing scattering, Phys. Rev. Lett. 17 (1966) 1158-1161.

[31] Z. Wang, Q. Guo, W. Hong, W. Hu, Modulational instability in nonlocal Kerr media with sine-oscillatory response, 2016. ARXIV 1608.01822v1[nlin.PS].

[32] Y.R. Shen, Principle of Nonlinear Optics, 1st edition, New York John Wiley and Sons, 1987.

[33] J.A. Hassen, C. Penland, Efficient Approximate Techniques for Integrating Stochastic Differential Equations, Monthly weather rev. 134 (2006) 3006-3014.

[34] T. Sauer, in, J.-C. Duan, W.K. Härdle, J.E. Gentle (Eds.), Handbook of Computational Finance, Springer, Heidelberg, Dordrecht, London, New York, 2012, pp. 529-550. Chapter 19.

[35] D.J. Higham, An Algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (2001) 525-546.

[36] A.J. Roberts, Modify the Improved Euler scheme to integrate stochastic differential equations, 2012. ARXIV 1210.0933 [math.NA].

[37] Y.E. Geints, A.A. Zemlyanov, Single and multiple filamentation of multi-terawatt CO₂-laser pulses in air: numerical simulations, J. OPT. SOC. AM. B 31 (2014) 788-797.

[38] S. Mitryukovskiy, Y. Liu, P. Ding, A. Houard, A. Mysyrowicz, Backward stimulated radiation from filaments in nitrogen gas and air pumped by circularly polarized 800 nm femtosecond laser pulses, Opt. Express 22 (2014) 12750-12759.

[39] J. Schwarz, J.C. Diels, Analytical solution for uv filaments, Phys. Rev. A 65 (2001) 013806-1-013806-10.

[40] Y. Shimoji, A.T. Fay, R.S.F. Chang, N. Djeu, Direct measurement of the nonlinear refractive index of air, J. Opt. Soc. Am. B 6 (1989) 1994-1998.

[41] C. Ting, H. Xu, C. Lin, 1.06-micron high power laser propagation in low-altitude atmosphere, in: Proceedings of the SPIE, 12413, 2023, pp. 124130I-1-124130I-13.

[42] P. Qi, W. Qian, L. Guo, J. Xue, N. Zhang, Y. Wang, Z. Zhang, Z. Zhang, L. Lin, C. Sun, et al., Sensing with Femtosecond Laser, Filamentation, Sensors 22 (2022) 7076.

[43] D. Zhu, C. Li, X. Sun, Y. Liu, Y. Zhang, H. Gao, The Effect of Air Turbulence on Vortex Beams in Nonlinear Propagation, Sensors 23 (2023) 1772.