Turbulence-resistant self-focusing vortex beams

Structured light, and in particular, light carrying orbital angular momentum (OAM) has been the subject of intense research since the seminal work by Allen and Woerdman in 1992 [1]. Several fields of fundamental and applied research, such as microscopy [2], spectroscopy [3], particle manipulation [4], quantum information [5], and optical communications [6], may benefit from light that carries OAM. A comprehensive review of the field, and a discussion of its future challenges can be found in references [7, 8], respectively. Amongst the different possibilities offered by structured light, vortex beams are the most promising candidates for the realisation of optical communication protocols, since they possess a virtually infinite alphabet that allows dense coding, multiplexing, and high-bit-rate communications [9, 10]. Vortex beams possess a twisted wave front, at the core of which sits a topologically-protected phase singularity described by the characteristic helical wavefront e^imϕ , where ϕ is the azimuthal angle around the beam propagation axis. The index $m\in \mathbb{Z}$ is known as the topological charge, which defines the amount of OAM carried by the beam, i.e., ℏm [7]. The unbounded nature of the index m translates to an infinite dimensional Hilbert space available to OAM states as opposed to the bounded two-dimensional Hilbert space associated to polarisation, which is ultimately what makes OAM enticing for optical communication purposes.

One of the main limiting factors that prevents OAM-carrying beams to be efficiently employed for free space optical communication, is the signal degradation due to atmospheric turbulence [11, 12]. For instance, one significant challenge that needs to be resolved is the modal cross-talk caused by turbulence, which can dramatically impair system performance [13, 14]. The dispersion of the OAM spectrum comes from the aberrant wavefront distorted by turbulence. The stronger the turbulence is, the more the spectrum will disperse [15]. Additionally, vortex beams experience significant divergence upon propagation—which is made worse by turbulence—making detection difficult due to limitations in detector size at the receiver plane [16]. Mitigating the effects of turbulence is an active area of study, and several schemes have been proposed to overcome these limitations, including adaptive-optics-enabled turbulence compensation [17], Laguerre–Gaussian mode sorters [18], and even deep learning [19]. The use of self-healing fields such as Airy or Bessel beams [20] seems tempting for this purpose. However, employing such beams in long range telecommunications is challenging, since the width of the beam is directly proportional to the desired propagation length, i.e. longer propagation requires a larger beam.

However, all the cases mentioned above consider the propagation of a fully coherent beam through the atmosphere. It is widely known that partially coherent fields feature superior resilience to intensity signal degradation upon propagation through turbulent media [21–25]. These results hint at the possibility of utilising partially coherent vortex beams to implement free space optical communication channels. Indeed, it has been revealed that one can reduce the turbulence-induced scintillation by decreasing the coherence of a vortex beam [26], which comes at a price of increased power loss over the link. Meanwhile, other studies have suggested that partially coherent vortex beams perform worse in turbulence than their coherent counterparts [27, 28]. These studies considered Schell-model type correlations (i.e. correlations that depend only on the distance between points), which is very similar to the effect of turbulence to begin with. Partially coherent electromagnetic fields may also feature locally varying degree of coherence, such as in the case of nonuniformly correlated fields [29], which feature self-focusing upon propagation. It has been shown that self-focusing beams possess better resilience against turbulence-induced noise than fully coherent fields [25, 30].

In the present work, we consider a recently introduced class of partially coherent self-focusing fields carrying OAM [31] and show how the combination of the innate resilience of self-focusing and information carrying capacity of OAM makes them a viable candidate for the realization of an efficient, free space optical communication link. In particular, we show that this kind of vortex beam exhibits self-focusing both in free space and in oceanic turbulence. In fact, it is found that the self-focusing property may be more pronounced in stronger turbulence and that the normalized OAM density shows also focusing behavior, compared to the coherent counterpart. This translates to a beam with several orders of magnitude higher peak intensity at the receiver, combined with an increased detection probability of the sent OAM state.

To start our analysis, let us consider a statistically stationary scalar source located in the plane z = 0, radiating a beam-like field towards the positive half-space z > 0. The spatial coherence properties of a 2D source at points ${\boldsymbol{\rho }}_{1}^{\prime }=\left({x}_{1}^{\prime },{y}_{1}^{\prime }\right)$ and ${\boldsymbol{\rho }}_{2}^{\prime }=\left({x}_{2}^{\prime },{y}_{2}^{\prime }\right)$, at angular frequency ω, can be described by the cross-spectral density (CSD) function [32]

$$\begin{equation}{W}_{0}({\boldsymbol{\rho }}_{1}^{\prime },{\boldsymbol{\rho }}_{2}^{\prime },\omega )=\left\langle {E}_{0}^{\ast }({\boldsymbol{\rho }}_{1}^{\prime },\omega ){E}_{0}({\boldsymbol{\rho }}_{2}^{\prime },\omega )\right\rangle ,\end{equation} \tag{ 1 }$$

where E₀( ρ ', ω) is the complex electric field, asterisk denotes the complex conjugate, and the angular brackets stand for ensemble averaging. We consider explicitly monochromatic fields, and thus dispersion effects are neglected. In what follows, the frequency dependence of all quantities of interest will be left implicit for brevity of notation.

It is well-known from standard coherence theory, that for a CSD to be physically realizable, it must be expressible in the form [33]

$$\begin{equation}{W}_{0}({\boldsymbol{\rho }}_{1}^{\prime },{\boldsymbol{\rho }}_{2}^{\prime })={\int }_{-\infty }^{\infty }p(v){H}_{0}^{\ast }({\boldsymbol{\rho }}_{1}^{\prime },v){H}_{0}({\boldsymbol{\rho }}_{2}^{\prime },v)\mathrm{d}v,\end{equation} \tag{ 2 }$$

where p(v) is a non-negative, ${\mathrm{L}}^{2}({\mathbb{R}}^{2})$-integrable probability density function, and H₀( ρ ', v) is a positive semidefinite kernel which we can choose as

$$\begin{equation}{H}_{0}({\boldsymbol{\rho }}^{\prime },v)=\tau ({\boldsymbol{\rho }}^{\prime })\mathrm{exp}\left(-2\pi \mathrm{i}v{\rho }^{\prime 2}\right),\end{equation} \tag{ 3 }$$

without loss of generality. Here τ( ρ ') is a (possibly complex valued) amplitude profile function and ${\rho }^{\prime }=\left\vert {\boldsymbol{\rho }}^{\prime }\right\vert$. The parameter v is an auxiliary variable which may or may not depend on the source properties, and which must vanish when the integration of equation (2) is carried out. Thus, it is a mathematical tool and with a dimension that depends on the form of the kernels. For example, v has the dimension of spatial frequency for a Fourier kernel.

The corresponding CSD can then be expressed as

$$\begin{equation}{W}_{0}(\boldsymbol{\rho }{;}_{1},{\boldsymbol{\rho }}_{2}^{\prime })={\tau }^{\ast }({\boldsymbol{\rho }}_{1}^{\prime })\tau ({\boldsymbol{\rho }}_{2}^{\prime }){\mu }_{0}({\rho }_{1}^{\prime 2}-{\rho }_{2}^{\prime 2}),\end{equation} \tag{ 4 }$$

where μ₀ is the (generalized) complex degree of spectral coherence, corresponding to the Fourier transform of p(v). Note that the profile function $\tau ({\boldsymbol{\rho }}_{1,2}^{\prime })$ is responsible only for the resulting intensity distribution, whereas the weight function p(v) produces the correlations. For the special case of v = 0, the above expression reduces to the CSD of a fully coherent field, which we will take as a reference, to benchmark the performance of the self-focusing OAM beams in the next sections. It is also worth mentioning, that since the only assumption on the weight function p(v) is to be square integrable, i.e., Fourier-transformable, regardless of the explicit form of p(v), equation (4) will always produce a self-focusing field.

From now on we make the explicit choice for the function p(v) to be a Gaussian,

$$\begin{equation}p(v)={\sigma }^{2}\sqrt{\pi }\,\mathrm{exp}\left(-{\pi }^{2}{\sigma }^{4}{v}^{2}\right)\end{equation} \tag{ 5 }$$

where σ indicates the coherence width of the beam. Beams with a larger value of σ have longer focal distances, which may be even on the order of kilometers [30]. Furthermore, we take the intensity profile of a Laguerre Gaussian beam with radial index zero [7, 34]

$$\begin{equation}\tau ({\boldsymbol{\rho }}^{\prime })={\left(\frac{{\rho }^{\prime }}{{w}_{0}}\right)}^{\vert m\vert }\,\mathrm{exp}(\mathrm{i}m\phi )\mathrm{exp}\left(-\frac{{\rho }^{\prime 2}}{2{w}_{0}^{2}}\right),\end{equation} \tag{ 6 }$$

where w₀ is the transverse beam width at z = 0, and m is the topological charge. By inserting this into equation (4), we end up with

$$\begin{align}\hfill {W}_{0}({\boldsymbol{\rho }}_{1}^{\prime },{\boldsymbol{\rho }}_{2}^{\prime })& ={\left(\frac{{\rho }_{1}^{\prime }{\rho }_{2}^{\prime }}{{w}_{0}^{2}}\right)}^{\vert m\vert }\,\mathrm{exp}[\mathrm{i}m({\phi }_{2}-{\phi }_{1})]\mathrm{exp}\left(-\frac{{\rho }_{1}^{\prime 2}+{\rho }_{2}^{\prime 2}}{2{w}_{0}^{2}}\right)\mathrm{exp}\left[-\frac{{({\rho }_{2}^{\prime 2}-{\rho }_{1}^{\prime 2})}^{2}}{{\sigma }^{4}}\right],\hfill \end{align} \tag{ 7 }$$

The negative sign in the last exponent, refers to the fact that we have a converging field. Therefore, equation (7) represents a self-focusing vortex beam, carrying m units of OAM.

Under the paraxial approximation, the propagation of a partially coherent beam from the source plane z = 0 to an arbitrary plane z > 0 in a turbulent medium can be described by the extended Huygens–Fresnel integral [32, 35]

$$\begin{align}\hfill W({\boldsymbol{\rho }}_{1},{\boldsymbol{\rho }}_{2},z)& ={\left(\frac{k}{2\pi z}\right)}^{2}{\iint }_{-\infty }^{\infty }{W}_{0}({\boldsymbol{\rho }}_{1}^{\prime },{\boldsymbol{\rho }}_{2}^{\prime })\mathrm{exp}\left[-\mathrm{i}k\frac{{\left({\boldsymbol{\rho }}_{1}-{\boldsymbol{\rho }}_{1}^{\prime }\right)}^{2}-{\left({\boldsymbol{\rho }}_{2}-{\boldsymbol{\rho }}_{2}^{\prime }\right)}^{2}}{2z}\right]\hfill \\ \hfill & \quad \times {\left\langle \mathrm{exp}\left[\varphi ({\boldsymbol{\rho }}_{1},{\boldsymbol{\rho }}_{1}^{\prime },z)+{\varphi }^{\ast }({\boldsymbol{\rho }}_{2},{\boldsymbol{\rho }}_{2}^{\prime },z)\right]\right\rangle }_{\text{M}}{\mathrm{d}}^{2}{\boldsymbol{\rho }}_{1}^{\prime }{\mathrm{d}}^{2}{\boldsymbol{\rho }}_{2}^{\prime },\hfill \end{align} \tag{ 8 }$$

where ${\boldsymbol{\rho }}_{1}=\left({x}_{1},{y}_{1}\right)$ and ${\boldsymbol{\rho }}_{2}=\left({x}_{2},{y}_{2}\right)$ represent two arbitrary spatial positions in the target plane, k = 2π/λ is the wave number, φ( ρ , ρ ', z) denotes the phase perturbation induced by the refractive-index fluctuations of the random medium between ρ ' and ρ , and ${\left\langle \dots \right\rangle }_{\text{M}}$ is the ensemble average over M realizations of turbulent media. This equation takes into account the random phase fluctuations, but neglects absorption. Since we are investigating monochromatic beams, the absorption can be added later by hand, if necessary.

If the fluctuations of the medium are homogeneous and isotropic, then the ensemble average in equation (8) can be expressed analytically as [35] (see also [36])

$$\begin{align}\hfill & {\langle \mathrm{exp}\left[{\varphi }^{\ast }({\boldsymbol{\rho }}_{1},{\boldsymbol{\rho }}_{1}^{\prime },z)+\varphi ({\boldsymbol{\rho }}_{2},{\boldsymbol{\rho }}_{2}^{\prime },z)\right]\rangle }_{\text{M}}\hfill \\ \hfill & \hspace{25.0pt}=\mathrm{exp}\left\{-\frac{1}{3}{\pi }^{2}{k}^{2}zT\left[{({\boldsymbol{\rho }}_{1}-{\boldsymbol{\rho }}_{2})}^{2}+({\boldsymbol{\rho }}_{1}-{\boldsymbol{\rho }}_{2})\cdot ({\boldsymbol{\rho }}_{1}^{\prime }-{\boldsymbol{\rho }}_{2}^{\prime })+{({\boldsymbol{\rho }}_{1}^{\prime }-{\boldsymbol{\rho }}_{2}^{\prime })}^{2}\right]\right\},\hfill \end{align} \tag{ 9 }$$

where T is the turbulence parameter. Here, we consider oceanic turbulence to find the performance of the beam in the worst conditions, in which case T has the explicit form [37]

$$\begin{equation}T=0.388\times 1{0}^{-8}{\varepsilon }^{-1/3}{\chi }_{\mathrm{T}}\left(47.5708{\varpi }^{-2}-17.6701{\varpi }^{-1}+6.783\,35\right).\end{equation} \tag{ 10 }$$

In the above equation, ɛ is the rate of dissipation of turbulent kinetic energy per unit mass of fluid, which may vary in the range [10⁻¹, 10⁻¹⁰ m² s⁻³], χ_T is the rate of dissipation of mean-square temperature, which has values in the range [10⁻⁴, 10⁻¹⁰] K² s⁻¹. Furthermore, the (dimensionless) parameter ϖ denotes the relative strength of temperature and salinity fluctuations, which for ocean water, ranges from −5 to 0, whose limits correspond to dominating temperature-induced and salinity-induced optical turbulence, respectively. The parameter T usually attains values ranging from about 10⁻¹⁶ m⁻¹ for weak turbulence, up to 10⁻¹² m⁻¹, for very strong turbulence although it has no well-defined upper limit [30]. Since the individual parameters in equation (10) are not relevant to our study, we use the value of T directly, mainly in the moderate to strong region of 10⁻¹⁴ m⁻¹ to 10⁻¹² m⁻¹.

To obtain the expression of the propagated CSD in presence of turbulence, we then substitute equation (2) into (8), such that at an arbitrary plane z we have the following CSD

$$\begin{equation}W({\boldsymbol{\rho }}_{1},{\boldsymbol{\rho }}_{2},z)={\int }_{-\infty }^{\infty }p(v){H}^{\ast }({\boldsymbol{\rho }}_{1},v,z)H({\boldsymbol{\rho }}_{2},v,z)\mathrm{d}v,\end{equation} \tag{ 11 }$$

where the propagated kernels are given by

$$\begin{align}\hfill {H}^{\ast }({\boldsymbol{\rho }}_{1},v,z)H({\boldsymbol{\rho }}_{2},v,z)& ={\left(\frac{k}{2\pi z}\right)}^{2}{\iint }_{-\infty }^{\infty }{H}_{0}^{\ast }({\boldsymbol{\rho }}_{1}^{\prime },v){H}_{0}({\boldsymbol{\rho }}_{2}^{\prime },v)\mathrm{exp}\left[-\mathrm{i}k\frac{{\left({\boldsymbol{\rho }}_{1}-{\boldsymbol{\rho }}_{1}^{\prime }\right)}^{2}-{\left({\boldsymbol{\rho }}_{2}-{\boldsymbol{\rho }}_{2}^{\prime }\right)}^{2}}{2z}\right]\hfill \\ \hfill & \quad \times {\langle \mathrm{exp}\left[{\varphi }^{\ast }({\boldsymbol{\rho }}_{1},{\boldsymbol{\rho }}_{1}^{\prime },z)+\varphi ({\boldsymbol{\rho }}_{2},{\boldsymbol{\rho }}_{2}^{\prime },z)\right]\rangle }_{\text{M}}{\mathrm{d}}^{2}{\boldsymbol{\rho }}_{1}^{\prime }{\mathrm{d}}^{2}{\boldsymbol{\rho }}_{2}^{\prime }.\hfill \end{align} \tag{ 12 }$$

Upon substituting equations (3) and (6) into equation (12) and going through some lengthy but straightforward algebra, we obtain the following expression for the product of two propagation kernels of a self-focusing vortex beam

$$\begin{align}\hfill {H}^{\ast }({\boldsymbol{\rho }}_{1},v,z)H({\boldsymbol{\rho }}_{2},v,z)& =\mathrm{exp}\left[-{\Theta}(x,y,T)-\frac{ik}{2z}\left({\boldsymbol{\rho }}_{1}^{2}-{\boldsymbol{\rho }}_{2}^{2}\right)-{T}_{1}{\left({\boldsymbol{\rho }}_{1}-{\boldsymbol{\rho }}_{2}\right)}^{2}\right]\hfill \\ \hfill & \quad \times {w}_{m}^{2}(z)\sum\limits _{n=0}^{\vert m\vert }{\left({\tilde{C}}_{\vert m\vert }^{n}\right)}^{2}{\left[F(x,y,T)\right]}^{\vert m\vert -n}.\hfill \end{align} \tag{ 13 }$$

Here T₁ = π² k² zT/3 is the z-dependent turbulence parameter, and we have employed multiple shorthand notations, all of which can be found in appendix A together with a detailed derivation of the equation. By employing the propagated kernels for each v together with equation (11), we can find the correlation properties of the self-focusing vortex beam at any plane of propagation.

4.1. Spectral density

We now use the formulae derived in the previous section to investigate the properties of the self-focusing vortex beams. When the two points coincide, ρ ₁ = ρ ₂ = ρ , equation (13) describes the spectral density of a single v-mode, and can be written in a simplified form as

$$\begin{equation}S(\boldsymbol{\rho },v,z)={w}_{m}^{2}(z)\mathrm{exp}\left[-\frac{{\boldsymbol{\rho }}^{2}}{{w}^{2}(z)}\right]\sum\limits _{n=0}^{\vert m\vert }{\left({\tilde{C}}_{\vert m\vert }^{n}\right)}^{2}{\left[\frac{k{w}_{\mathrm{f}\mathrm{r}}(z)\boldsymbol{\rho }}{zw(z)}\right]}^{2(\vert m\vert -n)}.\end{equation} \tag{ 14 }$$

For the case T = 0, i.e., in absence of turbulence, the waist of a single v-mode, w(z), reduces to its free space counterpart w_fr(z) (see appendix A). Further, if the OAM index is set to m = 0, then the intensity distribution reduces to a self-focusing field with no angular momentum, such as in reference [30], i.e., $S(\boldsymbol{\rho },v,z)={[{w}_{0}/w(z)]}^{2}\mathrm{exp}[-{\boldsymbol{\rho }}^{2}/w{(z)}^{2}]$. The overall spectral density [32] can be retrieved with

$$\begin{equation}S(\boldsymbol{\rho },z)=W(\boldsymbol{\rho },\boldsymbol{\rho },z)={\int }_{-\infty }^{\infty }p(v)S(\boldsymbol{\rho },v,z)\mathrm{d}v.\end{equation} \tag{ 15 }$$

In the following, we analytically simulate a beam propagating in oceanic turbulence. The employed parameters are given in the figure captions, whereas the mode waist, w₀ = 1 cm, and coherence width, σ = 5 mm are kept constant. These choices are not the only possible, and they produce a beam with about 60 m focal distance. Moreover, the wavelength was chosen to be 632.8 nm, which represents a sub-optimal choice. In a real-world application the operation wavelength should be chosen to minimize scattering and absorption, i.e. blue-green part of the visible spectrum when operating under water [38–40]. We chose a wavelength outside this window to test the model in worst-case conditions. As the sign of m just changes the handedness of the phase, we consider the case of positive topological charge only. Free space propagation can be regarded as the simplest case where the turbulence term as well as the summation index are both null, i.e. T = 0, n = 0 (see appendix A).

Figure 1 shows the evolution of the spectral density of the self-focusing vortex beam in the ρ − z plane in free space (a), and oceanic turbulence, for both weak (b), and strong (c) turbulence. Moreover, panel (d) shows the transverse distribution of the overall spectral density S( ρ , z) for different values of m at a selected distance z = 60 m. The peak value of the intensity and the radius of the dark zone increase for larger value of the topological charge, as expected for OAM-carrying beams [7]. A comparison between self-focusing vortex beam and a fully coherent vortex beam is made in panel (e), to help to evaluate the focusing effect. Here, we show the first main result of our work: the peak intensity of the self-focusing vortex beam is more than five orders of magnitude larger than that of a fully coherent vortex beam, although their intensities are of the same magnitude at the initial plane. This is due to the correlation induced self-focusing.

Zoom In Zoom Out Reset image size

4.2. Mean squared beam width

Let us next evaluate how the width of the beam evolves as it propagates in a turbulent ocean, by defining the mean squared beam width as

$$\begin{equation}{w}_{\mathrm{m}\mathrm{s}}^{2}(z)=\frac{{\int }_{0}^{\infty }{\boldsymbol{\rho }}^{2}S(\boldsymbol{\rho },z)\mathrm{d}\boldsymbol{\rho }}{{\int }_{0}^{\infty }S(\boldsymbol{\rho },z)\mathrm{d}\boldsymbol{\rho }}.\end{equation} \tag{ 16 }$$

Substituting equation (14) into the above expression, one can easily see that the behavior of ${w}_{\mathrm{m}\mathrm{s}}^{2}(z)$ is influenced by the OAM order m, the turbulence parameter T, and the propagation distance z, as one would expect. From this, one can determine the distance z_min, at which ${w}_{\mathrm{m}\mathrm{s}}^{2}$ reaches its minimum. However, since z_min involves complicated integrals and summations that do not necessarily admit closed-form solutions, we have to turn to numerical evaluation of the mean squared waist. We first consider a single component of the partially coherent beam, i.e., we take v as fixed. Using equations (14) and (16) we can evaluate the mean squared beam width for the vth component of the self-focusing vortex as

$$\begin{equation}{w}_{\mathrm{m}\mathrm{s}}^{2}(v,z)={w}^{2}(v,z){Y}_{m}\left[1-\frac{{w}^{2}(v,z)}{{w}_{\mathrm{f}\mathrm{r}}^{2}(v,z)}\right],\end{equation} \tag{ 17 }$$

where the v dependence has been written explicitly for clarity, and

$$\begin{equation}{Y}_{m}(x)=\frac{{\Gamma}\left(\frac{3}{2}+\vert m\vert \right){\Gamma}{\left(\frac{1}{2}-\vert m\vert \right)}_{2}{F}_{1}\left(-\vert m\vert ,-\vert m\vert ,-\frac{1}{2}-\vert m\vert ,x\right)}{{\pi {(-1)}^{\vert m\vert }}_{2}{F}_{1}\left(-\vert m\vert ,-\vert m\vert ,\frac{1}{2}-\vert m\vert ,x\right)}\end{equation} \tag{ 18 }$$

is the OAM dependent correction term, with Γ(n) being the Gamma function, and ₂ F₁(a, b; c; x) the Gauss hypergeometric function [41].

Let us briefly mention some of the main features of the mean squared beam width: (i) for m = 0, we have Y₀(x) = 1 for any x. Moreover, in free space for any given |m| > 0 we have Y_m (0) = 1 as well. (ii) Since w(z) increases monotonically with T, the mean squared width of the vth component will also monotonically increase with turbulence when Y_m (x) = 1. (iii) The value of the function Y_m (−|x|) increases with increasing m, and thus, high OAM content beams have a larger mean square width. (iv) For a fixed |m| > 0, Y_m (x) decreases monotonically with the increase of turbulence strength, until it tends to a constant value after a sufficient propagation distance. These findings are graphically summarized in figure 2.

Zoom In Zoom Out Reset image size

It is evident that the mean squared beam width becomes narrower over a certain range of distances, as the turbulence increases (dashed lines in figure 2). This becomes more obvious for beams with larger OAM content, such as those depicted in panels (c) and (d). Additionally, one can see how the influence of the term Y_m (x) in equation (17) affects the focusing dynamics of the beam by modifying the minimum value, as well as the focal distance.

Previous studies of partially coherent self-focusing beams (without vortices) in turbulent ocean has revealed that increasing turbulence causes the focal spot to move towards the source plane [30]. However, in the case of self-focusing vortex beams, the picture is somewhat different. In this case, the focal spot moves towards to the source plane only for low values of OAM, as it can be seen in figures 2(a) and (b) for m = 1. For higher values of OAM, the position of z_min is shifted away from the source plane, as shown in figures 2(c) and (d) for m = 8 and m = 12.

This constitutes the second main result of our work: by evaluating the mean squared beam width we find that, for higher OAM content and within certain distances, the spread of a partially coherent vortex beam in strong turbulence can be lower than in weak turbulence. This is, to the best of our knowledge, the first example of an optical beam possessing such a counter-intuitive property.

4.3. Orbital angular momentum flux density

We now proceed to our analysis of the OAM properties of the beam. For a scalar and paraxial partially coherent beam, which has on average no spin angular momentum, the z component of the OAM flux density is of the form [42]

$$\begin{equation}{L}_{\text{orb}}(\boldsymbol{\rho },z)=-\frac{{{\epsilon}}_{0}}{k}\mathrm{I}\mathrm{m}{\left[{y}_{1}\frac{\partial }{\partial {x}_{2}}W({\boldsymbol{\rho }}_{1},{\boldsymbol{\rho }}_{2},z)-{x}_{1}\frac{\partial }{\partial {y}_{2}}W({\boldsymbol{\rho }}_{1},{\boldsymbol{\rho }}_{2},z)\right]}_{{\boldsymbol{\rho }}_{1}={\boldsymbol{\rho }}_{2}},\end{equation} \tag{ 19 }$$

where ₀ denotes the free space permittivity. The quantity L_orb is dependent on the vortex strength, as well as the intensity of the field.

To have a better understanding of the physical picture, we consider the normalized OAM flux density, which describes the OAM per photon [42],

$$\begin{equation}{l}_{\text{orb}}(\boldsymbol{\rho },z)=\frac{\hslash \omega {L}_{\text{orb}}(\boldsymbol{\rho },z)}{{S}_{\mathrm{p}}(\boldsymbol{\rho },z)},\end{equation} \tag{ 20 }$$

where S_p( ρ , z) = kW( ρ , ρ , z)/(μ₀ ω) is the z component of the Poynting vector, with μ₀ being the vacuum permeability. Upon substituting equations (11) and (13) into equation (20), the normalized OAM flux density can be written as

$$\begin{equation}{l}_{\text{orb}}(\boldsymbol{\rho },z)=\hslash \frac{{\int }_{-\infty }^{\infty }\,p\,(v){S}_{l}(\boldsymbol{\rho },v,z)\mathrm{d}v}{{\int }_{-\infty }^{\infty }\,p\,(v)S(\boldsymbol{\rho },v,z)\mathrm{d}v},\end{equation} \tag{ 21 }$$

where

$$\begin{align}\hfill {S}_{l}(\boldsymbol{\rho },v,z)& =\frac{{w}_{m}^{2}(z){w}^{2}(z)}{{w}_{\mathrm{f}\mathrm{r}}^{2}(z)}\mathrm{exp}\left[-\frac{{\boldsymbol{\rho }}^{2}}{{w}^{2}(z)}\right]\sum\limits _{n=0}^{\vert m\vert }{({\tilde{C}}_{\vert m\vert }^{n})}^{2}(\vert m\vert -n){\left[\frac{k{w}_{\mathrm{f}\mathrm{r}}(z)\boldsymbol{\rho }}{zw(z)}\right]}^{2(\vert m\vert -n)},\hfill \end{align} \tag{ 22 }$$

Using equations (21) and (22), we can study the behaviour of the normalized OAM flux density as a function of the radial position at any propagation distance of a self-focusing vortex beam carrying m units of OAM.

By comparing the expressions in equations (14) and (22), one can readily see that in free space (i.e., T = 0) equation (21) reduces to the well known constant value of l_orb( ρ , z) = |m|ℏ. Thus, self-focusing vortex beams in free space can be viewed as pure fluid-like rotators for which the OAM content is independent of the radial distance. Therefore, the type of self-focusing vortex beams we consider act like pure coherent vortex modes in free space with respect to their OAM content [43, 44], since l_orb( ρ , z) is simply proportional to the topological charge m.

When turbulence is present, the OAM flux density attains a nontrivial dependence on radial, as well as axial distance. The effect of turbulence on the OAM flux density of partially coherent self-focusing vortex beams for different m at selected distances are shown in figure 3.

Zoom In Zoom Out Reset image size

Compared to the coherent vortex beams, the self-focusing vortex beams behave more like Rankine vortices, with a small rigid body rotation region near the core (where the OAM depends quadratically on the distance ρ) and a fluid-like rotation in the outer regions. With increasing turbulence and propagation distance, the rigid region expands. This is not surprising since, as it follows directly from the last term in equation (22), ρ is inversely proportional to z and T. Intriguingly, the maximum value of the OAM flux density for the self-focusing vortex beam exceeds the value of the topological charge carried by the input beam, and it takes a higher value under stronger turbulence. This can be explained as the result of balancing between turbulence induced spreading and the self-focusing of the beam.

4.4. Orbital angular momentum spectrum

Finally, we investigate the effects of oceanic turbulence on the OAM spectrum of self-focusing vortex beams. It is well known that there is a Fourier relationship between the intensity distribution in the azimuthal direction and the complex OAM spectrum [45], i.e.,

$$\begin{equation}E(\rho ,\theta ,z)=\frac{1}{\sqrt{2\pi }}\sum\limits _{l=-\infty }^{+\infty }\,{a}_{l}(\rho ,z)\mathrm{exp}(\mathrm{i}l\theta ),\end{equation} \tag{ 23 }$$

and

$$\begin{equation}{a}_{l}(\rho ,z)=\frac{1}{\sqrt{2\pi }}{\int }_{0}^{2\pi }E(\rho ,\theta ,z)\mathrm{exp}(-\mathrm{i}l\theta )\mathrm{d}\theta .\end{equation} \tag{ 24 }$$

The power weights of partially coherent vortex beams can be written as

$$\begin{equation}\left\langle \vert {a}_{l}(\rho ,z){\vert }^{2}\right\rangle =\frac{1}{2\pi }{\iint }_{0}^{2\pi }W(\rho ,{\theta }_{1},\rho ,{\theta }_{2},z)\mathrm{exp}[\mathrm{i}l({\theta }_{1}-{\theta }_{2})]\mathrm{d}{\theta }_{1}\mathrm{d}{\theta }_{2},\end{equation} \tag{ 25 }$$

in which case the relative power for each OAM mode is defined as

$$\begin{equation}{P}_{l}(z)=\frac{{\int }_{0}^{\infty }\left\langle \vert {a}_{l}(\rho ,z){\vert }^{2}\right\rangle \rho \mathrm{d}\rho }{{\sum }_{l}{\int }_{0}^{\infty }\left\langle \vert {a}_{l}(\rho ,z){\vert }^{2}\right\rangle \rho \mathrm{d}\rho }.\end{equation} \tag{ 26 }$$

Computing the relative powers of the states with various azimuthal indices allows us to construct the OAM spectrum [46].

The OAM spectrum of a self-focusing vortex beam, as well as its coherent counterpart for various propagation distances and turbulence parameters are displayed in figure 4. Since we saw from the previous subsection that the OAM flux density is a constant in free space for these beams, we can immediately surmise that the OAM spectrum in free space contains only a single peak at m (the signal mode). Therefore, in the following we consider only beams in turbulence.

Zoom In Zoom Out Reset image size

With increasing propagation distance and turbulence parameter, the energy of the signal mode will inject into the neighbouring modes for both self-focusing and coherent vortex beams, resulting in OAM spectrum dispersion. Thus, the probability of detecting the signal mode decreases monotonously upon propagation through turbulence. The partially coherent self-focusing beam has a heavier tailed OAM spectrum than the fully coherent one. However, the proportion of the signal channel is higher than that of coherent ones, featuring a maximum increase of ∼14% over the coherent vortex beam. This is partly caused by the narrower beam size of the self-focusing beam, which ensures smaller area affected by turbulence.

As a last example of the present work, we study the OAM spectrum received by a detector at two propagation distances, which is depicted in figure 5. Here, we assume that the receiver plane may move along the z-axis, and that the employed detector is circular, with a radius of 10 mm. Further, we assume perfect alignment, and thus the detector area contains most of the energy of the self-focusing beam.

Zoom In Zoom Out Reset image size

This brings us to the third main result of our work. By comparing figures 4 and 5, we see that the finite aperture cuts the tails of the OAM spectrum. Hence, the weight of the signal channel becomes much larger, especially for the case of moderately strong oceanic turbulence $( > 0.8)$. Furthermore, the proportion of the overall energy that the detector captures from the self-focusing vortex beam with m = 2 in moderately strong turbulence is 72.3% at z = 60 m and 67.7% at z = 150 m. In the case of strong turbulence, these values are 72.9% and 68.9%, respectively. Comparing this to a coherent vortex beam, the ratios at the same propagation distances are 15% and 14.7% with T = 10⁻¹⁴ m⁻¹, whereas they are 15.1% and 16.5% with T = 10⁻¹³ m⁻¹, respectively. The energy efficiency in the last coherent example is higher in stronger turbulence due to an increase of intensity on axis (i.e. due to vortex breakdown). In each considered case, there is an over four-fold increase in efficiency when compared to a completely coherent beam, and the improvement is greater for longer distances when the focal distance is modified accordingly.

We can also see that the probability of the signal mode decreases with the increase of topological charge m, indicating that a vortex beam with a greater m is more vulnerable to turbulence, as larger m gives rise to a wider beam radius, which consequently causes serious OAM dispersion. This means that although it is a good idea to construct high dimensional Hilbert space by increasing the value of the topological charge m, it may not be efficient in improving free space optical communication. In addition, it is to be noted that the OAM spectrum becomes asymmetric and shifts towards zero under increasing turbulence, which is due to the beam losing angular momentum to the turbulence [47]. Mathematically, this can be explained from the turbulence correlated coordinates a_s and ξ_s, as the shifts of the coordinates in x and y directions are asymmetric (see appendix A for more information).

We have studied the propagation properties of self-focusing vortex beams through oceanic turbulence. We chose ocean environment to find the worst case estimates on the propagation properties of such beams, and the results are applicable also to the much weaker atmospheric turbulence. The self-focusing vortex beams feature several intriguing properties, and they can be extremely robust against turbulence.

Although the present work is entirely theoretical, it is of interest to discuss how such beams could be experimentally generated. Circularly correlated beams—which also feature self-focusing—have been realized with digital micro-mirror devices [48], where an input beam was modulated with a time dependent hologram. Nonuniform correlations (the correlation type we study) may be generated in a similar fashion, or with deformable optics. In the latter case, the optical element would be a mirror with variable curvature, driven with a suitable time dependent signal. After the correlations have been engineered, one only needs to impose an OAM phase on the output beam to attain the self-focusing vortex beam. This may be done with standard holographic techniques, which can also generate high-order OAM beams. Moreover, the detection of these partially coherent beams can be done in the same fashion as coherent beams. In particular, for detecting higher-order OAM beams, one could again use standard holographic techniques based on SLMs, or OAM sorter techniques (see, for example, references [49, 50]) to detect the complete OAM spectrum of the beam after propagation through turbulence. Alternatively, one may also measure the recently introduced coherence OAM matrix [51], which can be done with a modified Young's interferometer [52].

First, thanks to the focusing of the beam upon propagation, the peak intensity at the receiver plane is over five orders of magnitude greater than for a corresponding coherent beam. This is an important feature for a beam employed in a data link, since it makes alignment and detection much easier. Here, we considered a beam which has a focal length of about 60 m. However, by increasing the coherence width of the beam, it is possible to get much longer focal distances, even on the order of kilometers. To retain the robustness of the beam, the overall degree of coherence has to remain unchanged, and thus an increase of the coherence width necessitates an increase in the overall beam size. Since the size of the beam is finite in real world applications, this limits the longest achievable focal length. For example, with beam and coherence widths of 20 mm, one can already get a focal length of 0.8 km [30].

Second, the width of the beam and the focal distance remain nearly unaltered when the turbulence is increased. In fact, under certain conditions, the mean squared beam waist may become smaller for a larger turbulence parameter. This effect appears to be due to an interplay between OAM and turbulence mediated by the focusing properties of the beam. We believe that we are the first to predict such counter-intuitive behaviour in optical beams.

Third, the signal mode detection probability of the partially coherent beam is improved over the coherent one, and when evaluating OAM spectrum over the whole beam, the signal mode detection probability increases by about 14%. When a circular detector of 10 mm radius is employed at the receiving plane, the signal mode is detected with a probability of about 80% over moderately strong turbulence. Moreover, the finite detector can collect ∼70% of the energy that is sent over a distance of ∼100 m through strong turbulence, making for an extremely efficient link.

In addition to the properties listed above, partial coherence leads to a significant decrease in scintillation, and therefore improves the receiver's signal-to-noise ratio. Thus, a data link realized with self-focusing OAM beams has several advantages over coherent OAM beams. The results obtained in the present study may find applications in long distance OAM data links, whether it is through an atmospheric or oceanic environment.

ML acknowledges the financial support of the National Natural Science Foundation of China (NSFC) (61805080), Hunan Provincial Natural Science Foundation of China (2019JJ50366) and the China Scholarship Council (CSC). MK, and MO acknowledge the financial support of the Academy of Finland Flagship Programme (PREIN-decisions 320165, 320166). CD acknowledges the financial support of the National Natural Science Foundation of China (NSFC) (12174171) and 2020 Central Plains Talents Program of Henan (ZYYCYU202012144).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In the following, we assume that the initial plane correlation function is of the form of a self-focusing vortex beam that is given in equation (7), the derivation of which was outlined in section 2 of the main text. Further, we assume paraxial propagation towards the positive half-space z > 0, which is described by the extended Huygens–Fresnel integral of equation (8). By combining these, we attain the kernel propagation integral of equation (12). Upon substituting from equations (3) and (6) into equation (12), we get

$$\begin{align}\hfill {H}^{\ast }({\boldsymbol{\rho }}_{1},v,z)H({\boldsymbol{\rho }}_{2},v,z)& ={\left(\frac{k}{2\pi z}\right)}^{2}\iint {\left({\rho }_{1}^{\prime }{\rho }_{2}^{\prime }\right)}^{\vert m\vert }{\mathrm{e}}^{-im({\phi }_{1}-{\phi }_{2})}\mathrm{exp}\left(-\frac{{\rho }_{1}^{\prime 2}+{\rho }_{2}^{\prime 2}}{2{w}_{0}^{2}}\right)\mathrm{exp}\left[-2\pi iv\left({\rho }_{2}^{\prime 2}-{\rho }_{1}^{\prime 2}\right)\right]\hfill \\ \hfill & \quad \times \mathrm{exp}\left\{-\frac{ik}{2z}\left[{\left({\boldsymbol{\rho }}_{1}-{\boldsymbol{\rho }}_{1}^{\prime }\right)}^{2}-{\left({\boldsymbol{\rho }}_{2}-{\boldsymbol{\rho }}_{2}^{\prime }\right)}^{2}\right]\right.\hfill \\ \hfill & \quad -{T}_{1}\left[{\left({\boldsymbol{\rho }}_{1}-{\boldsymbol{\rho }}_{2}\right)}^{2}+({\boldsymbol{\rho }}_{1}-{\boldsymbol{\rho }}_{2})\cdot ({\boldsymbol{\rho }}_{1}^{\prime }-{\boldsymbol{\rho }}_{2}^{\prime })\left.+{({\boldsymbol{\rho }}_{1}^{\prime }-{\boldsymbol{\rho }}_{2}^{\prime })}^{2}\right]\right\}{\mathrm{d}}^{2}{\boldsymbol{\rho }}_{1}^{\prime }{\mathrm{d}}^{2}{\boldsymbol{\rho }}_{2}^{\prime },\hfill \end{align} \tag{ A1 }$$

where T₁ = π² k² zT/3 is the z-dependent turbulence parameter.

To make the integral simpler, we can separate it to ${\boldsymbol{\rho }}_{1}^{\prime }$ and ${\boldsymbol{\rho }}_{2}^{\prime }$ contributions. The integral over ${\boldsymbol{\rho }}_{1}^{\prime }$ takes the following—analytically solvable—form in Cartesian coordinates

$$\begin{align}\hfill & \iint {\left[{x}_{1}^{\prime }-i\mathrm{sgn}(m){y}_{1}^{\prime }\right]}^{\vert m\vert }\mathrm{exp}\left[-{\alpha }_{1}\left({x}_{1}^{\prime 2}+{y}_{1}^{\prime 2}\right)+{u}_{x}{x}_{1}^{\prime }+{u}_{y}{y}_{1}^{\prime }\right]\mathrm{d}{x}_{1}^{\prime }\mathrm{d}{y}_{1}^{\prime }\hfill \\ \hfill & \hspace{25.0pt}=\frac{\pi }{{\alpha }_{1}}{\left[\frac{{u}_{x}-i\mathrm{sgn}(m){u}_{y}}{2{\alpha }_{1}}\right]}^{\vert m\vert }\mathrm{exp}\left(\frac{{u}_{x}^{2}+{u}_{y}^{2}}{4{\alpha }_{1}}\right),\hfill \end{align} \tag{ A2 }$$

where sgn(x) is the sign function. Here we have used the shorthand notations

$$\begin{equation}{\alpha }_{1}=\frac{1}{2{w}_{0}^{2}}-i2\pi v+\frac{ik}{2z}+{T}_{1},\end{equation} \tag{ A3 }$$

$$\begin{equation}{u}_{x}=\frac{ik{x}_{1}}{z}+{T}_{1}{\Delta}x+2{T}_{1}{x}_{2}^{\prime },\end{equation} \tag{ A4 }$$

$$\begin{equation}{u}_{y}=\frac{ik{y}_{1}}{z}+{T}_{1}{\Delta}y+2{T}_{1}{y}_{2}^{\prime },\end{equation} \tag{ A5 }$$

$$\begin{equation}{\Delta}x={x}_{2}-{x}_{1},\hspace{25.0pt}{\Delta}y={y}_{2}-{y}_{1}.\end{equation} \tag{ A6 }$$

Next, we choose a_x = kx₁/z − iT₁Δx and a_y = ky₁/z − iT₁Δy, so that we can write ${u}_{x}=i{a}_{x}+2{T}_{1}{x}_{2}^{\prime }$, ${u}_{y}=i{a}_{y}+2{T}_{1}{y}_{2}^{\prime }$, then the remaining integral over ${\boldsymbol{\rho }}_{2}^{\prime }$ takes the form

$$\begin{align}\hfill & \iint {[i{a}_{x}+\mathrm{sgn}(m){a}_{y}+2{T}_{1}({x}_{2}^{\prime }-i\mathrm{sgn}(m){y}_{2}^{\prime })]}^{\vert m\vert }{({x}_{2}^{\prime }+i\mathrm{sgn}(m){y}_{2}^{\prime })}^{\vert m\vert }\mathrm{exp}\left[-\left({\alpha }_{1}^{\ast }-\frac{{T}_{1}^{2}}{{\alpha }_{1}}\right)\right.\hfill \\ \hfill & \hspace{25.0pt}\times \left({x}_{2}^{\prime 2}+{y}_{2}^{\prime 2}\right)-\left(\frac{ik}{z}{x}_{2}+{T}_{1}{\Delta}x-\frac{i{a}_{x}{T}_{1}}{{\alpha }_{1}}\right){x}_{2}^{\prime }\left.-\left(\frac{ik}{z}{y}_{2}+{T}_{1}{\Delta}y-\frac{i{a}_{y}{T}_{1}}{{\alpha }_{1}}\right){y}_{2}^{\prime }\right]\mathrm{d}{x}_{2}^{\prime }\mathrm{d}{y}_{2}^{\prime },\hfill \end{align} \tag{ A7 }$$

Using the formula

$$\begin{equation}{(a+b)}^{m}=\sum\limits _{n=0}^{m}\,{C}_{m}^{n}{a}^{m-n}{b}^{n},\end{equation} \tag{ A8 }$$

where ${C}_{m}^{n}=m!/n!(m-n)!$, we can expand ${[i{a}_{x}+\mathrm{sgn}(m){a}_{y}+2{T}_{1}({x}_{2}^{\prime }-i\mathrm{sgn}(m){y}_{2}^{\prime })]}^{\vert m\vert }$ in a binomial series to get

$$\begin{align}\hfill & \sum\limits _{n=0}^{\vert m\vert }{C}_{\vert m\vert }^{n}{[i{a}_{x}+\mathrm{sgn}(m){a}_{y}]}^{(\vert m\vert -n)}{(2{T}_{1})}^{n}\iint {\left[{x}_{2}^{\prime }-i\mathrm{sgn}(m){y}_{2}^{\prime }\right]}^{n}{\left[{x}_{2}^{\prime }+i\mathrm{sgn}(m){y}_{2}^{\prime }\right]}^{\vert m\vert }\hfill \\ \hfill & \hspace{25.0pt}\times \mathrm{exp}\left[-\left({\alpha }_{1}^{\ast }-\frac{{T}_{1}^{2}}{{\alpha }_{1}}\right)\left({x}_{2}^{\prime 2}+{y}_{2}^{\prime 2}\right)-i{\xi }_{x}{x}_{2}^{\prime }-i{\xi }_{y}{y}_{2}^{\prime }\right]\mathrm{d}{x}_{2}^{\prime }\mathrm{d}{y}_{2}^{\prime },\hfill \end{align} \tag{ A9 }$$

where we have used

$$\begin{equation}{\xi }_{x}=\frac{k{x}_{2}}{z}-i{T}_{1}{\Delta}x-\frac{{a}_{x}{T}_{1}}{{\alpha }_{1}},\end{equation} \tag{ A10 }$$

$$\begin{equation}{\xi }_{y}=\frac{k{y}_{2}}{z}-i{T}_{1}{\Delta}y-\frac{{a}_{y}{T}_{1}}{{\alpha }_{1}}.\end{equation} \tag{ A11 }$$

Now, we can write equation (A9) in cylindrical coordinates as

$$\begin{equation}\begin{aligned}\hfill & \sum\limits _{n=0}^{\vert m\vert }{C}_{\vert m\vert }^{n}{[i{a}_{x}+\text{sgn}\;(m){a}_{y}]}^{(\vert m\vert -n)}{(2{T}_{1})}^{n}\int {{\rho }_{2}^{\prime }}^{\vert m\vert +n+1}\mathrm{exp}\left[-\left({\alpha }_{1}^{\ast }-\frac{{T}_{1}^{2}}{{\alpha }_{1}}\right){\rho }_{2}^{\prime 2}\right]\mathrm{d}{\rho }_{2}^{\prime }\hfill \\ \hfill & \quad \times \int \mathrm{exp}\left[\mathrm{sgn}(m)i(\vert m\vert -n){\phi }_{2}^{\prime }-i{\xi }_{x}{\rho }_{2}^{\prime }\cos {\phi }_{2}^{\prime }-i{\xi }_{y}{\rho }_{2}^{\prime }\mathrm{sin}{\phi }_{2}^{\prime }\right]\mathrm{d}{\phi }_{2}^{\prime }\hfill \end{aligned}\end{equation} \tag{ A12 }$$

and further employ the formula [53]

$$\begin{equation}{\int }_{0}^{2\pi }{\mathrm{e}}^{\pm \mathrm{i}\mathrm{n}\mathrm{t}+i(x\mathrm{sin}t+y\cos t)}\mathrm{d}t=\frac{2\pi {J}_{\pm n}\left(\sqrt{{x}^{2}+{y}^{2}}\right)}{{(-1)}^{n}{\sqrt{{x}^{2}+{y}^{2}}}^{n}}{(x\mp iy)}^{n},\end{equation} \tag{ A13 }$$

to solve the azimuthal angle integral, as well as [53]

$$\begin{equation}{\int }_{0}^{\infty }{x}^{2n+u+1}{\mathrm{e}}^{-{x}^{2}}{J}_{u}(2x\sqrt{z})\mathrm{d}x=\frac{n!}{2}{\mathrm{e}}^{-z}{z}^{\frac{1}{2}u}{L}_{n}^{u}(z),\end{equation} \tag{ A14 }$$

$$\begin{equation}{L}_{n}^{\alpha }(x)=\sum\limits _{k=0}^{n}\frac{{(-1)}^{k}}{k!}{C}_{n+\alpha }^{n-k}{x}^{k},\end{equation} \tag{ A15 }$$

for the radial integral. After some calculations, the product of the propagated kernels is obtained as

$$\begin{align}\hfill {H}^{\ast }({\boldsymbol{\rho }}_{1},v,z)H({\boldsymbol{\rho }}_{2},v,z)& =\mathrm{exp}\left[-{\Theta}(x,y,T)-\frac{ik}{2z}\left({\boldsymbol{\rho }}_{1}^{2}-{\boldsymbol{\rho }}_{2}^{2}\right)-{T}_{1}{\left({\boldsymbol{\rho }}_{1}-{\boldsymbol{\rho }}_{2}\right)}^{2}\right]\hfill \\ \hfill & \quad \times {w}_{m}^{2}(z)\sum\limits _{n=0}^{\vert m\vert }{\left({\tilde{C}}_{\vert m\vert }^{n}\right)}^{2}{\left[F(x,y,T)\right]}^{\vert m\vert -n}.\hfill \end{align} \tag{ A16 }$$

which is equation (13) in the main text, and the employed shorthand notations are

$$\begin{equation}{\Theta}(x,y,T)=({a}_{x}^{2}+{a}_{y}^{2})/(4{\alpha }_{1})+({\xi }_{x}^{2}+{\xi }_{y}^{2})/(4{\alpha }_{2}),\end{equation} \tag{ A17 }$$

$$\begin{equation}{w}_{m}(z)={z}^{\vert m\vert }{w}_{0}/[{k}^{\vert m\vert }w{(z)}^{\vert m\vert +1}],\end{equation} \tag{ A18 }$$

$$\begin{equation}{\tilde{C}}_{\vert m\vert }^{n}=\frac{\vert m\vert !}{n!(\vert m\vert -n)!}\sqrt{n!{(4{T}_{1})}^{n}},\end{equation} \tag{ A19 }$$

$$\begin{equation}F(x,y,T)=[{a}_{y}+\mathrm{sgn}(m)i{a}_{x}][{\xi }_{y}-\mathrm{sgn}(m)i{\xi }_{x}]-({T}_{1}/{\alpha }_{2})({\xi }_{x}^{2}+{\xi }_{y}^{2}),\end{equation} \tag{ A20 }$$

$$\begin{equation}{\alpha }_{2}={\alpha }_{1}^{\ast }-\frac{{T}_{1}^{2}}{{\alpha }_{1}}.\end{equation} \tag{ A21 }$$

Lastly, the turbulence-corrected, z-dependent waist of the beam is given by

$$\begin{equation}{w}^{2}(z)={\left(\frac{z}{k{w}_{0}}\right)}^{2}+4{T}_{1}{\left(\frac{z}{k}\right)}^{2}+{\left(1-\frac{4\pi vz}{k}\right)}^{2}{w}_{0}^{2},\end{equation} \tag{ A22 }$$

whereas the free space waist is simply

$$\begin{equation}{w}_{\mathrm{f}\mathrm{r}}^{2}(z)={\left(\frac{z}{k{w}_{0}}\right)}^{2}+{\left(1-\frac{4\pi vz}{k}\right)}^{2}{w}_{0}^{2}.\end{equation} \tag{ A23 }$$

Employing these equations together with the main text considerations, allows us to analytically model a self-focusing vortex beam of arbitrary OAM content at any propagation distance.