The behavior of partially coherent twisted space-time beams in atmospheric turbulence

We study how atmospheric turbulence affects twisted space-time beams, which are non-stationary random optical fields whose space and time dimensions are coupled with a stochastic twist. Applying the extended Huygens–Fresnel principle, we derive the mutual coherence function of a twisted space-time beam after propagating a distance z through atmospheric turbulence of arbitrary strength. We specialize the result to derive the ensemble-averaged irradiance and discuss how turbulence affects the beam’s spatial size, pulse width, and space-time twist. Lastly, we generate, in simulation, twisted space-time beam field realizations and propagate them through atmospheric phase screens to validate our analysis.

1 Introduction

New approaches in beam control include light with engineered space-time or spatiotemporal coupling. Recent papers have demonstrated space-time-coupled light which exhibits anomalous diffractive and refractive behaviors [1–4] as well as carries transverse (to the direction of propagation) orbital angular momentum in the form of spatiotemporal optical vortices (STOVs) [5–11]. These novel developments hold promise for exciting advancements in applications such as optical communications, optical tweezing, and quantum optics [2, 4, 12–16].

Most of the space-time-coupled beam research manipulates coherent light, although this has begun to change with the development of partially coherent STOV and twisted space-time (and space-frequency) beams [17–21]. The latter are non-stationary random fields with the beams’ spatial and temporal (or spectral) dimensions coupled in a stochastic twist. They are the spatiotemporal counterparts of traditional, spatially twisted Gaussian Schell-model beams [22–27].

Spatially twisted partially coherent fields have been extensively studied since being introduced in 1993. This research includes beam synthesis [28–33]; coherent modes/pseudo-modes [23, 26, 27, 34–37]; angular momentum [38–41]; and propagation behaviors in free-space, ABCD optical systems, and turbulence [35, 42–51]. This stands in contrast to twisted space-time beams (and STOV beams more generally), where only their angular momentum and free-space propagation behaviors have been investigated [6, 8, 9, 11, 19, 20, 52, 53].

In this paper, we undertake, to our knowledge, the first study on the effects of atmospheric turbulence on twisted space-time beams. Using the extended Huygens–Fresnel principle, we derive an approximate expression for the mutual coherence function (MCF) of a twisted space-time beam after propagating through atmospheric turbulence of any strength. We then specialize the MCF to obtain the ensemble-averaged irradiance and discuss how turbulence affects the beam’s size, pulse width, and space-time twist. To validate our analysis, we compare the theoretical irradiance to the results of Monte Carlo wave-optics simulations. Lastly, we conclude with a brief summary of our findings.

2 Theory

2.1 Extended Huygens–Fresnel principle

Let us begin with the extended Huygens–Fresnel principle/integral:

$U\left(\mathit{\rho },z,\omega \right)=\frac{k}{\text{j}2\pi z}\exp \left(\text{j}kz\right){\iint }_{-\infty }^{\infty }U\left({\mathit{\rho }}^{\prime },0,\omega \right)\exp \left(\frac{\text{j}k}{2z}{\left|\mathit{\rho }-{\mathit{\rho }}^{\prime }\right|}^2\right)\exp \left[\mathrm{\Psi }\left({\mathit{\rho }}^{\prime },0;\mathit{\rho },z;\omega \right)\right]{\text{d}}^2{\rho }^{\prime },(1)$

where $\text{j}=\sqrt{-1}$, ω is the radian frequency, k = ω/c is the wavenumber, c is the speed of light, ${\mathit{\rho }}^{\prime }=\widehat{\mathit{x}}{x}^{\prime }+\widehat{\mathit{y}}{y}^{\prime }$ is the source vector, and $\mathit{\rho }=\widehat{\mathit{x}}x+\widehat{\mathit{y}}y$ is the observation vector. The optical field U in the integrand is a stochastic (frequency-domain) realization of a twisted space-time beam, and Ψ is a random complex function which models the phase and amplitude fluctuations of a point source propagating through atmospheric turbulence from $\left({\mathit{\rho }}^{\prime },0\right)$ to $\left(\mathit{\rho },z\right)$ at frequency ω [54–57].

The two-frequency cross-spectral density (CSD) function [55, 58–62] can be obtained by taking the ensemble-averaged auto-correlation of Eq. 1, namely,

$\begin{array}{c}W\left({\mathit{\rho }}_1,z,{\omega }_1,{\mathit{\rho }}_2,z,{\omega }_2\right)=\frac{k_1{k}_2}{{\left(2\pi \right)}^2{z}^2}\exp \left[\text{j}\left(k_1-k_2\right)z\right]{\iint \iint }_{-\infty }^{\infty }W\left({\mathit{\rho }}_1^{\prime },0,{\omega }_1,{\mathit{\rho }}_2^{\prime },0,{\omega }_2\right)\\ \exp \left(\frac{\text{j}{k}_1}{2z}{\left|{\mathit{\rho }}_1-{\mathit{\rho }}_1^{\prime }\right|}^2\right)\exp \left(-\frac{\text{j}{k}_2}{2z}{\left|{\mathit{\rho }}_2-{\mathit{\rho }}_2^{\prime }\right|}^2\right)\\ ⟨\exp \left[\mathrm{\Psi }\left({\mathit{\rho }}_1^{\prime },0;{\mathit{\rho }}_1,z;{\omega }_1\right)+{\mathrm{\Psi }}^{*}\left({\mathit{\rho }}_2^{\prime },0;{\mathit{\rho }}_2,z;{\omega }_2\right)\right]⟩{\text{d}}^2{\rho }_1^{\prime }{\text{d}}^2{\rho }_2^{\prime },\end{array}(2)$

where we have assumed that the source field is statistically independent of the atmospheric turbulence fluctuations. The moment involving Ψ is related to the two-point, spherical wave structure function (WSF) [55–57, 61, 62], and equals

$\begin{array}{c}⟨\exp \left[\mathrm{\Psi }\left({\mathit{\rho }}_1^{\prime },0;{\mathit{\rho }}_1,z;{\omega }_1\right)+{\mathrm{\Psi }}^{*}\left({\mathit{\rho }}_2^{\prime },0;{\mathit{\rho }}_2,z;{\omega }_2\right)\right]⟩=\exp \left[-\frac{1}{2}D\left({\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime },0;{\mathit{\rho }}_1-{\mathit{\rho }}_2,z;{\omega }_1,{\omega }_2\right)\right]\\ \approx \exp \left\{-2{\pi }^2{\int }_0^z{\int }_0^{\infty }\kappa {\mathrm{\Phi }}_n\left(\kappa ,\zeta \right)\left[k_1^2+k_2^2-2k_1{k}_2\exp \left(-\text{j}\beta {\kappa }^2\right)J_0\left(\kappa R\right)\right]\text{d}\kappa \text{d}\zeta \right\},\end{array}(3)$

where Φ_n is the index of refraction power spectrum (assumed to be statistically isotropic) and β and R equal

$\begin{array}{cc}\hfill \beta & =\frac{\zeta \left(z-\zeta \right)}{2z}\left(\frac{1}{k_1}-\frac{1}{k_2}\right)\hfill \\ \hfill R& =\left|\frac{\zeta }{z}\left({\mathit{\rho }}_1-{\mathit{\rho }}_2\right)+\left(1-\frac{\zeta }{z}\right)\left({\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime }\right)\right|.\hfill \end{array}(4)$

The approximate expression on the second line of Eq. 3 is derived using the method of smooth perturbations (also known as the Rytov approximation) and further assuming that Ψ is Gaussian distributed [55–57, 61–63]. We will return to Eq. 3 shortly.

The ultimate goal is to find the “two-time” MCF of a twisted space-time beam after propagating through turbulence. To do this, we must inverse Fourier transform Eq. 2, i.e.,

$\mathrm{\Gamma }\left({\mathit{\rho }}_1,z,t_1,{\mathit{\rho }}_2,z,t_2\right)={\iint }_{-\infty }^{\infty }W\left({\mathit{\rho }}_1,z,{\omega }_1,{\mathit{\rho }}_2,z,{\omega }_2\right)\exp \left(-\text{j}{\omega }_1{t}_1\right)\exp \left(\text{j}{\omega }_2{t}_2\right)\text{d}{\omega }_1\text{d}{\omega }_2.(5)$

Applying Eqs 2–5 and interchanging the order of the integrals yields

$\begin{array}{c}\mathrm{\Gamma }\left({\mathit{\rho }}_1,z,t_1,{\mathit{\rho }}_2,z,t_2\right)\\ =\frac{1}{{\left(2\pi \right)}^2{c}^2{z}^2}{\iint \iint }_{-\infty }^{\infty }{\iint }_{-\infty }^{\infty }{\omega }_1{\omega }_{2}W\left({\mathit{\rho }}_1^{\prime },0,{\omega }_1,{\mathit{\rho }}_2^{\prime },0,{\omega }_2\right)\\ \exp \left[-\frac{1}{2}D\left({\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime },0;{\mathit{\rho }}_1-{\mathit{\rho }}_2,z;{\omega }_1,{\omega }_2\right)\right]\\ \exp \left[-\text{j}{\omega }_1\left(t_1-\frac{z}{c}-\frac{{\left|{\mathit{\rho }}_1-{\mathit{\rho }}_1^{\prime }\right|}^2}{2zc}\right)\right]\\ \exp \left[\text{j}{\omega }_2\left(t_2-\frac{z}{c}-\frac{{\left|{\mathit{\rho }}_2-{\mathit{\rho }}_2^{\prime }\right|}^2}{2zc}\right)\right]\text{d}{\omega }_1\text{d}{\omega }_2{\text{d}}^2{\rho }_1^{\prime }{\text{d}}^2{\rho }_2^{\prime }.\end{array}(6)$

Assuming that the twisted space-time beam has a relatively narrow linewidth (or bandwidth) around mean or carrier frequency ω_c (i.e., Δω/ω_c ≪ 1), we can approximate Eq. 6 as

$\begin{array}{c}\mathrm{\Gamma }\left({\mathit{\rho }}_1,z,t_1,{\mathit{\rho }}_2,z,t_2\right)\approx \frac{1}{{\lambda }_c^2{z}^2}{\iint \iint }_{-\infty }^{\infty }\exp \left[-\text{j}{\omega }_c\left(t_1-t_2-\frac{{\left|{\mathit{\rho }}_1-{\mathit{\rho }}_1^{\prime }\right|}^2}{2zc}+\frac{{\left|{\mathit{\rho }}_2-{\mathit{\rho }}_2^{\prime }\right|}^2}{2zc}\right)\right]\\ {\iint }_{-\infty }^{\infty }W\left({\mathit{\rho }}_1^{\prime },0,{\bar{\omega }}_1,{\mathit{\rho }}_2^{\prime },0,{\bar{\omega }}_2\right)\exp \left[-\frac{1}{2}D\left({\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime },0;{\mathit{\rho }}_1-{\mathit{\rho }}_2,z;{\bar{\omega }}_1+{\omega }_c,{\bar{\omega }}_2+{\omega }_c\right)\right]\exp \left[-\text{j}{\bar{\omega }}_1\left(t_1-\frac{z}{c}\right)\right]\exp \left[\text{j}{\bar{\omega }}_2\left(t_2-\frac{z}{c}\right)\right]\text{d}{\bar{\omega }}_1\text{d}{\bar{\omega }}_2{\text{d}}^2{\rho }_1^{\prime }{\text{d}}^2{\rho }_2^{\prime },\end{array}(7)$

and, by evaluating Eq. 7, obtain a closed-form expression for the MCF. Before doing this, we need to discuss the functions D and W in the integrand.

2.2 Approximate two-point, spherical WSF D

Let us return to Eq. 3. By virtue of the source being narrowband, β ≈ 0. Letting Φ_n equal the von Kármán spectrum—namely,

${\mathrm{\Phi }}_n\left(\kappa ,\zeta \right)=0.033C_n^2\left(\zeta \right)\frac{\exp \left(-{\kappa }^2/{\kappa }_m^2\right)}{{\left(\kappa +{\kappa }_0^2\right)}^{11/6}},(8)$

where κ_m = 5.92/l₀ and κ₀ = 2π/L₀ ($C_n^2$, l₀, and L₀ are the index of refraction structure constant, the inner scale, and outer scale of turbulence, respectively)—the integral over κ evaluates to

$\begin{array}{c}D\left({\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime },0;{\mathit{\rho }}_1-{\mathit{\rho }}_2,z;{\omega }_1,{\omega }_2\right)=0.033\left(2{\pi }^2\right){\kappa }_0^{-5/3}U\left(1;\frac{1}{6};\frac{{\kappa }_0^2}{{\kappa }_m^2}\right)\\ \times \left(k_1^2+k_2^2\right){\int }_0^z{C}_n^2\left(\zeta \right)\text{d}\zeta \\ -0.033\left(4{\pi }^2\right){\kappa }_0^{-5/3}{k}_1{k}_2{\int }_0^z{C}_n^2\left(\zeta \right)\\ \times \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{n!}{\left(\frac{{\kappa }_0^2{R}^2}{4}\right)}^{n}U\left(n+1;n+\frac{1}{6};\frac{{\kappa }_0^2}{{\kappa }_m^2}\right)\text{d}\zeta ,\end{array}(9)$

where $U\left(a;c;z\right)$ is a confluent hypergeometric function of the second kind [64–66]. In most physical scenarios, L₀ ≫ l₀, and therefore, we can estimate Eq. 9 using the small argument relation for $U\left(a;c;z\right)$. The result, after much analysis, is

$\begin{array}{c}D\left({\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime },0;{\mathit{\rho }}_1-{\mathit{\rho }}_2,z;{\omega }_1,{\omega }_2\right)\approx 0.7817{\kappa }_0^{-5/3}{\left(k_1-k_2\right)}^2{\int }_0^z{C}_n^2\left(\zeta \right)\text{d}\zeta \\ +8.7021{\kappa }_m^{-5/3}{k}_1{k}_2{\int }_0^z{C}_n^2\left(\zeta \right)\left[{}_{1}F_1\left(-\frac{5}{6};1;-\frac{{\kappa }_m^2{R}^2}{4}\right)-1\right]\text{d}\zeta -2.3450{\kappa }_0^{1/3}{k}_1{k}_2{\int }_0^z{C}_n^2\left(\zeta \right)R^2\text{d}\zeta .\end{array}(10)$

Eq. 10 includes both inner and outer scale effects. However, to evaluate Eq. 7 in closed form, we must let the inner scale l₀ → 0 (κ_m → ∞). Using the large argument relation for ₁F₁, we obtain

$\begin{array}{c}D\left({\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime },0;{\mathit{\rho }}_1-{\mathit{\rho }}_2,z;{\omega }_1,{\omega }_2\right)\approx 0.7817{\kappa }_0^{-5/3}{\left(k_1-k_2\right)}^2{\int }_0^z{C}_n^2\left(\zeta \right)\text{d}\zeta \\ +2.9139k_1{k}_2{\int }_0^z{C}_n^2\left(\zeta \right)R^{5/3}\text{d}\zeta -2.3450{\kappa }_0^{1/3}{k}_1{k}_2{\int }_0^z{C}_n^2\left(\zeta \right)R^2\text{d}\zeta .\end{array}(11)$

We lastly assume that $C_n^2$ is constant over the propagation path and set R^5/3 ≈ R²—an estimate known as the quadratic approximation [57, 62]. Substituting ${\omega }_1={\bar{\omega }}_1+{\omega }_c$ and ${\omega }_2={\bar{\omega }}_2+{\omega }_c$ as stipulated in Eq. 7 and noting that $k_1{k}_2=\left({\bar{k}}_1+k_c\right)\left({\bar{k}}_2+k_c\right)\approx k_c^2$, we arrive at the final result

$\begin{array}{c}D\left({\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime },0;{\mathit{\rho }}_1-{\mathit{\rho }}_2,z;{\bar{\omega }}_1+{\omega }_c,{\bar{\omega }}_2+{\omega }_c\right)\\ \approx \frac{0.7817C_n^{2}z{\kappa }_0^{-5/3}}{c^2}{\left({\bar{\omega }}_1-{\bar{\omega }}_2\right)}^2+1.0930C_n^{2}z\left(1-0.7152{\kappa }_0^{1/3}\right)k_c^2\\ \left[{\left|{\mathit{\rho }}_1-{\mathit{\rho }}_2\right|}^2+{\left|{\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime }\right|}^2+\left({\mathit{\rho }}_1-{\mathit{\rho }}_2\right)\cdot \left({\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime }\right)\right]\\ =2\frac{a_{\omega }}{c^2}{\left({\bar{\omega }}_1-{\bar{\omega }}_2\right)}^2+2a_s{k}_c^2\left[{\left|{\mathit{\rho }}_1-{\mathit{\rho }}_2\right|}^2+{\left|{\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime }\right|}^2+\left({\mathit{\rho }}_1-{\mathit{\rho }}_2\right)\cdot \left({\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime }\right)\right].\end{array}(12)$

Equation 12 is very physical: The terms describe how atmospheric turbulence corrupts light’s spectral and spatial coherence. For traditional space-time separable beams, these two terms give rise to pulse and beam broadening, respectively [56, 57, 67–73]. In our case, because of spatiotemporal coupling, both terms will affect the temporal and spatial beam sizes.

2.3 CSD function of a twisted space-time beam

With Eq. 12, we are one step closer to evaluating Eq. 7. We, of course, still need an expression for W. To find this expression, we begin with the MCF of a twisted space-time beam:

$\begin{array}{c}\mathrm{\Gamma }\left({\mathit{\rho }}_1,t_1,{\mathit{\rho }}_2,t_2\right)=A^2\exp \left(-\frac{y_1^2+y_2^2}{4W_y^2}\right)\exp \left(-\frac{x_1^2+x_2^2}{4W_x^2}\right)\exp \left[-\frac{{\left(x_1-x_2\right)}^2}{2{\delta }_x^2}\right]\\ \exp \left(-\frac{t_1^2+t_2^2}{4W_t^2}\right)\exp \left[-\frac{{\left(t_1-t_2\right)}^2}{2{\delta }_t^2}\right]\exp \left[\text{j}\mu \left(x_1{t}_2-x_2{t}_1\right)\right]\exp \left[-\text{j}{\omega }_c\left(t_1-t_2\right)\right];\end{array}(13)$

where A is the amplitude; W_x, W_y, and W_t are the spatial and temporal pulse widths; δ_x and δ_t are the spatial and temporal coherence widths; and μ is the space-time twist parameter [19]. The latter must satisfy $\left|\mu \right|{\delta }_t{\delta }_x\le 1$ for the MCF in Eq. 13 to be genuine, i.e., square-integrable, Hermitian, and non-negative definite [58, 59]. Consequently, μ → 0 in the coherent beam limit δ_t, δ_x → ∞. When $\left|\mu \right|{\delta }_t{\delta }_x=1$, the twist in the beam is saturated [20, 25, 27]. We assume this condition for the simulations described in Section 3.

Note that Eq. 13 has the same general form as a twisted Gaussian Schell-model beam [22, 25–27]; however, here, space and time are twisted. It is well known that the spectral density or average irradiance of a spatially twisted random beam rotates in the x-y plane as it propagates in the z direction [35, 40, 41, 74]. From Eq. 7, we see that t is linked paraxially to the propagation distance z; therefore, a twisted space-time beam rotates or tumbles in the x-z plane as it propagates. This behavior is described in Refs. [19, 20] for twisted space-time beams propagating in free space. What remains to be determined is how atmospheric turbulence affects the x-z plane rotation of twisted space-time beams.

We can find the CSD function W of a twisted space-time beam by Fourier transforming the MCF in Eq. 13, i.e.,

$W\left({\mathit{\rho }}_1,{\omega }_1,{\mathit{\rho }}_2,{\omega }_2\right)=\frac{1}{{\left(2\pi \right)}^2}{\iint }_{-\infty }^{\infty }\mathrm{\Gamma }\left({\mathit{\rho }}_1,t_1,{\mathit{\rho }}_2,t_2\right)\exp \left(\text{j}{\omega }_1{t}_1\right)\exp \left(-\text{j}{\omega }_2{t}_2\right)\text{d}{t}_1\text{d}{t}_2.(14)$

Substituting Eq. 13 into Eq. 14 and evaluating the integrals yields

$\begin{array}{c}W\left({\mathit{\rho }}_1,{\omega }_1,{\mathit{\rho }}_2,{\omega }_2\right)=\frac{A^2}{4\pi {\mathrm{\Omega }}_t}\exp \left(-\frac{y_1^2+y_2^2}{4W_y^2}\right)\exp \left(-\frac{x_1^2+x_2^2}{4{\tilde{W}}_x^2}\right)\exp \left[-\frac{{\left(x_1-x_2\right)}^2}{2{\tilde{\delta }}_x^2}\right]\\ \exp \left(-\frac{{\bar{\omega }}_1^2+{\bar{\omega }}_2^2}{4W_{\omega }^2}\right)\exp \left[-\frac{{\left({\bar{\omega }}_1-{\bar{\omega }}_2\right)}^2}{2{\delta }_{\omega }^2}\right]\exp \left\{\frac{\mu }{{\delta }_{\omega }^2}\left[{\bar{\omega }}_2\left(1+\frac{1}{2{\gamma }_t^2}\right)-{\bar{\omega }}_1\right]x_1\right\}\\ \exp \left\{\frac{\mu }{{\delta }_{\omega }^2}\left[{\bar{\omega }}_1\left(1+\frac{1}{2{\gamma }_t^2}\right)-{\bar{\omega }}_2\right]x_2\right\},\end{array}(15)$

where γ_t = W_t/δ_t, W_ω = 2W_tΩ_t, δ_ω = 2δ_tΩ_t, and

$\begin{array}{c}{\mathrm{\Omega }}_t^2={\left(\frac{1}{4W_t^2}+\frac{1}{2{\delta }_t^2}\right)}^2-{\left(\frac{1}{2{\delta }_t^2}\right)}^2,\\ \frac{1}{{\tilde{W}}_x^2}=\frac{1}{W_x^2}+\frac{{\mu }^2}{W_{\omega }^2},\\ \frac{1}{{\tilde{\delta }}_x^2}=\frac{1}{{\delta }_x^2}+\frac{{\mu }^2}{{\delta }_{\omega }^2}.\end{array}(16)$

With Eq. 15, we are now ready to evaluate the integrals in Eq. 7.

2.4 MCF of twisted space-time beam in atmospheric turbulence

Substituting Eqs 12, 15 into Eq. 7 and evaluating the integrals produces (after much analysis)

$\begin{array}{c}\mathrm{\Gamma }\left(x_1,0,z,t_1,x_2,0,z,t_2\right)=A^2\frac{W_t}{{\tilde{W}}_t}\frac{N_{Fx}{N}_{Fy}}{\sqrt{{\mathrm{\Delta }}_x{\mathrm{\Delta }}_y}}\exp \left[-\text{j}{\omega }_c\left({\bar{t}}_1-{\bar{t}}_2\right)\right]\\ \exp \left(-\frac{N_{Fx}^2}{{\mathrm{\Delta }}_x}\frac{x_1^2+x_2^2}{4W_x^2}\right)\exp \left[-\frac{{\bar{t}}_1^2+{\bar{t}}_2^2}{4{\left(W_t^{\text{eff}}\right)}^2}\right]\exp \left[\mu \frac{N_{Fx}}{{\mathrm{\Delta }}_x}\frac{W_t^2}{{\tilde{W}}_t^2}\left(x_1{\bar{t}}_2+x_2{\bar{t}}_1\right)\right]\\ \exp \left[\frac{\text{j}{k}_c}{2z}\left(1-\frac{N_{Fx}^2}{{\mathrm{\Delta }}_x}+a_s{k}_c^2\frac{4W_x^2}{{\mathrm{\Delta }}_x}\right)\left(x_1^2-x_2^2\right)\right]\exp \left[\text{j}{\mu }^2\frac{N_{Fx}^2}{{\mathrm{\Delta }}_x}\frac{W_t^2}{{\tilde{W}}_t^2}\left({\bar{t}}_1^2-{\bar{t}}_2^2\right)\right]\\ \exp \left[\text{j}\mu \frac{N_{Fx}^2}{{\mathrm{\Delta }}_x}\frac{W_t^2}{{\tilde{W}}_t^2}\left(x_1{\bar{t}}_2-x_2{\bar{t}}_1\right)\right]\exp \left[-\text{j}\mu \frac{N_{Fx}^2}{{\mathrm{\Delta }}_x}\frac{a_{\omega }}{c^2{\tilde{W}}_t^2}\left(x_1+x_2\right)\left({\bar{t}}_1-{\bar{t}}_2\right)\right]\\ \exp \left[-\text{j}\mu \frac{2W_x^2}{{\mathrm{\Delta }}_x}\frac{W_t^2}{{\tilde{W}}_t^2}{a}_s{k}_c^2\left(x_1-x_2\right)\left({\bar{t}}_1+{\bar{t}}_2\right)\right]\\ \exp \left[-\frac{{\left(x_1-x_2\right)}^2}{2{\left({\delta }_x^{\text{eff}}\right)}^2}\right]\exp \left[-\frac{{\left({\bar{t}}_1-{\bar{t}}_2\right)}^2}{2{\left({\delta }_t^{\text{eff}}\right)}^2}\right]\exp \left[-\frac{\left(x_1-x_2\right)\left({\bar{t}}_1-{\bar{t}}_2\right)}{2{\left({\delta }_{xt}^{\text{eff}}\right)}^2}\right].\end{array}(17)$

Since the beam’s interesting behaviors occur in the x-t or x-z plane (the x and t dimensions are coupled), here, we present the MCF evaluated at y₁ = y₂ = 0. The undefined symbols in Eq. 17 are $N_{Fx,y}=2k_c{W}_{x,y}^2/z$, which are the x and y Fresnel numbers for a fully coherent Gaussian beam; $\bar{t}=t-z/c$ is the retarded time; ${\tilde{W}}_t^2=W_t^2+2a_{\omega }/c^2$; and

$\begin{array}{c}{\mathrm{\Delta }}_x=1+4\frac{W_x^2}{{\delta }_x^2}+N_{Fx}^2+8W_x^2\left(a_s{k}_c^2+{\mu }^2\frac{a_{\omega }}{c^2}\frac{W_t^2}{{\tilde{W}}_t^2}\right);\\ {\mathrm{\Delta }}_y=1+N_{Fy}^2+8W_y^2{a}_s{k}_c^2;\\ \\ \frac{1}{4{\left(W_t^{\text{eff}}\right)}^2}=\frac{1}{4{\tilde{W}}_t^2}+{\mu }^2\frac{W_x^2}{{\mathrm{\Delta }}_x}\frac{W_t^4}{{\tilde{W}}_t^4};\\ \frac{1}{2{\left({\delta }_t^{\text{eff}}\right)}^2}=\frac{1}{2{\delta }_t^2}+{\mu }^2\frac{W_x^2}{2{\mathrm{\Delta }}_x}\left({\mathrm{\Delta }}_x-N_{Fx}^2-\frac{W_t^4}{{\tilde{W}}_t^4}\right)+\frac{a_{\omega }}{c^2}\frac{1}{4W_t^2{\tilde{W}}_t^2};\\ \\ \frac{1}{2{\left({\delta }_x^{\text{eff}}\right)}^2}=\frac{N_{Fx}^2}{{\mathrm{\Delta }}_x}\frac{1}{2{\delta }_x^2}+a_s{k}_c^2\left(1+2\frac{N_{Fx}^2}{{\mathrm{\Delta }}_x}\right)-{\left(a_s{k}_c^2\right)}^2\frac{2W_x^2}{{\mathrm{\Delta }}_x^2}+{\mu }^2\frac{a_{\omega }}{c^2}\frac{N_{Fx}^2}{{\mathrm{\Delta }}_x}\frac{W_t^2}{{\tilde{W}}_t^2};\\ \frac{1}{2{\left({\delta }_{xt}^{\text{eff}}\right)}^2}=\mu \frac{N_{Fx}}{2{\mathrm{\Delta }}_x}\left({\mathrm{\Delta }}_x-N_{Fx}^2-\frac{W_t^2}{{\tilde{W}}_t^2}+a_s{k}_c^{2}4W_x^2\right).\end{array}(18)$

Eq. 17 is organized so that the terms can be physically interpreted: Starting at the top and ignoring the carrier $\exp \left[-\text{j}{\omega }_c\left({\bar{t}}_1-{\bar{t}}_2\right)\right]$, the amplitude term plus the first three exponentials comprise the ensemble-averaged irradiance (discussed in more detail below). The next (complex) exponentials on line 4 are the spatial and temporal chirps. These are followed by the space-time twist on lines 5 and 6. Lastly, the exponentials on line 7 model spatial and temporal coherence.

2.5 Average irradiance and physical discussion

The ensemble-averaged irradiance is found by evaluating Eq. 17 at equal space-time points, i.e.,

$\begin{array}{c}⟨I\left(x,0,z,t\right)⟩=\mathrm{\Gamma }\left(x,0,z,t,x,0,z,t\right)\\ =A^2\frac{W_t}{{\tilde{W}}_t}\frac{N_{Fx}{N}_{Fy}}{\sqrt{{\mathrm{\Delta }}_x{\mathrm{\Delta }}_y}}\exp \left(-\frac{N_{Fx}^2}{{\mathrm{\Delta }}_x}\frac{x^2}{2W_x^2}\right)\exp \left[-\left(1+{\mu }^2\frac{4W_x^2{W}_t^4}{{\mathrm{\Delta }}_x{\tilde{W}}_t^2}\right)\frac{{\bar{t}}^2}{2{\tilde{W}}_t^2}\right]\exp \left(\mu \frac{2N_{Fx}}{{\mathrm{\Delta }}_x}\frac{W_t^2}{{\tilde{W}}_t^2}x\bar{t}\right)\\ ={\widehat{A}}^2\exp \left(-\frac{x^2}{2{\widehat{W}}_x^2}\right)\exp \left(-\frac{{\bar{t}}^2}{2{\widehat{W}}_t^2}\right)\exp \left(\widehat{\mu }x\bar{t}\right).\end{array}(19)$

In order, the exponentials are the spatial beam shape, temporal beam (pulse) shape, and x-t plane rotation. The behavior of the beam can be understood by examining ${\widehat{W}}_x$, ${\widehat{W}}_t$, and $\widehat{\mu }$ versus Fresnel number and turbulence strength. Figure 1 shows these curves: (A) plots ${\widehat{W}}_x/W_x$, ${\widehat{W}}_t/W_t$, and $\widehat{\mu }/\mu$ over Fresnel numbers ranging from 100 (near field) to 0.01 (far field). The solid, dashed, dashed-dotted, and dotted traces show how these quantities evolve in free space ($C_n^2=0{\text{\hspace{0.17em}m}}^{-2/3}$) and atmospheric turbulence ($C_n^2=10^{-14}{\text{\hspace{0.17em}m}}^{-2/3}$ with L₀ = 10 m, 50 m, and 100 m), respectively. For the latter, the (weak turbulence) spherical wave scintillation indices [57], i.e.,

${\sigma }_I^2=0.5C_n^2{k}_c^{7/6}{z}^{11/6},(20)$

are annotated on the plot (centered on their corresponding Fresnel number) to show the strength of turbulence at that N_Fx. Figure 1B displays a zoomed-in view of ${\widehat{W}}_x/W_x$, ${\widehat{W}}_t/W_t$, and $\widehat{\mu }/\mu$ over the boxed region in (A), viz., 15 ≥ N_Fx ≥ 0.5. Lastly, the results depicted in Figure 1 apply to a twisted space-time beam with λ_c = 1 μm, W_x = 2 cm, δ_x = 0.9W_x, W_t = 1 ps, δ_t = 0.9W_t, and $\mu =1/\left({\delta }_x{\delta }_t\right)$.

FIGURE 1

FIGURE 1. (A) ${\widehat{W}}_x/W_x$, ${\widehat{W}}_t/W_t$, and $\widehat{\mu }/\mu$ from Eq. 19 versus Fresnel number N_Fx. (B) Zoomed-in view of ${\widehat{W}}_x/W_x$, ${\widehat{W}}_t/W_t$, and $\widehat{\mu }/\mu$ over the boxed region in (A). The solid, dashed, dashed-dotted, and dotted traces show the results in free space $(C_n^2=0{\text{\hspace{0.17em}m}}^{-2/3})$ and atmospheric turbulence ($C_n^2=10^{-14}{\text{\hspace{0.17em}m}}^{-2/3}$ with L₀=10 m, 50 m, and 100 m), respectively. The text annotations report the weak turbulence spherical wave scintillation index values ${\sigma }_I^2$ at the corresponding N_Fx. These results apply to a twisted space-time beam with parameter values equal to λ_c =1 μm, W_x =2 cm, δ_x =0.9W_x, W_t =1 ps, δ_t =0.9W_t, and $\mu =1/\left({\delta }_x{\delta }_t\right)$.

Starting with the free-space (solid) curves in Figure 1, we see that for N_Fx > 10, the twisted space-time beam is effectively in the source plane, with ${\widehat{W}}_x\approx W_x$, ${\widehat{W}}_t\approx W_t$, and $\widehat{\mu }<\mu /5$. Things begin to change for 10 > N_Fx > 1: Most noticeably, the beam grows significantly larger due to diffraction. Indeed, over this range, the beam expands nearly three times its original size in the x direction. In addition to ${\widehat{W}}_x$, the pulse width also changes in this region because of spatiotemporal coupling. Beginning around N_Fx ≈ 10, ${\widehat{W}}_t$ starts to contract (shorten) and continues to do so until N_Fx ≈ 1. This shortening of ${\widehat{W}}_t$ is met by an increase in $\widehat{\mu }$. When considered together, the result is a beam that rotates in the x-t (or x-z) plane—the beam effectively “trades” ${\widehat{W}}_t$ to do so. Lastly, for N_Fx < 1, ${\widehat{W}}_x$ continues to grow larger due to diffraction, ${\widehat{W}}_t$ asymptotes (the pulse width stops contracting), and $\widehat{\mu }$ falls rapidly toward zero. Physically, the twisted space-time beam is in the far zone, diffraction dominates, and the beam no longer rotates.

Examining the turbulence (dashed, dashed-dotted, and dotted) curves, we generally observe the same behavior; however, the beam’s evolution described above is effectively pushed to the left, i.e., toward higher Fresnel numbers. Where the separation between free-space (diffractive) and turbulence-induced behavior occurs (in other words, at what N_Fx), of course, depends on $C_n^2$ and L₀. Nevertheless, some general trends are evident and independent of turbulence strength:

1. The beam’s size ${\widehat{W}}_x$ asymptotically expands much more rapidly in turbulence than in free space (z³ vice z²) [69, 71–73].

2. After initially contracting, the pulse width ${\widehat{W}}_t$ lengthens and continues to grow longer. While this can clearly be seen in Figure 1, more insight can be gained by examining the mathematical expression for ${\widehat{W}}_t$, namely,

$\frac{1}{{\widehat{W}}_t}=\frac{1}{{\tilde{W}}_t}\left(1+{\mu }^2\frac{4W_x^2{W}_t^4}{{\mathrm{\Delta }}_x{\tilde{W}}_t^2}\right).(21)$

In atmospheric turbulence, the term containing the twist parameter μ tends to zero like z⁻² (in free space, the term asymptotes to a constant value). For large z, the result is therefore ${\widehat{W}}_t\sim {\tilde{W}}_t=W_t+2a_{\omega }/c^2$. The turbulence contribution to the pulse width grows linearly with z [57, 67, 68], thus explaining the increasing pulse width.

3. The x-t plane rotation $\widehat{\mu }$ decays much more rapidly in turbulence than in free space. Examining the mathematical relation for $\widehat{\mu }$ reveals that it approaches zero like $\widehat{\mu }\sim z^{-3}$ in turbulence (vice $\widehat{\mu }\sim z^{-1}$ in free space) as z → ∞.

3 Validation

In this section, we validate Eq. 19 by generating, in simulation, twisted space-time beam field realizations and propagating those realizations through atmospheric turbulence phase screens. Before presenting and analyzing the results, we discuss the simulation setup.

3.1 Simulation setup

3.1.1 Numbers of grid points, spacings, trials, etc.

In these wave-optics simulations, we generated and propagated twisted space-time beam field realizations through independent instances of atmospheric turbulence. The Fresnel numbers for these simulations were N_Fx = 10, 5, 2.5, and 1. For each N_Fx, we computed the ensemble-averaged irradiance $⟨I\left(x,0,z,t\right)⟩$ from 1,000 independent field and turbulence realizations. The source and observation planes were discretized using three-dimensional grids that were N_y × N_x × N_t = 1, 200 × 1, 200 × 128 with spacings equal to Δ_src = 1.58 mm, Δ_obs = 2.5 mm, and Δt = 0.0781 ps.

3.1.2 Generating twisted space-time fields

We generated twisted space-time beam field realizations using the approach described in Ref. [31]. The technique utilizes Gori and Santarsiero’s integral criterion for genuine CSD functions and MCFs, colloquially known as the superposition rule [75, 76]. Specialized for our purposes, a thermal (or pseudo-thermal) twisted space-time beam field realization can be generated by evaluating the following superposition integral:

$U\left(\mathit{\rho },t\right)={\iint }_{-\infty }^{\infty }r\left(v_x,v_t\right)\sqrt{\frac{1}{2}p\left(v_x,v_t\right)}H\left(\mathit{\rho },t;v_x,v_t\right)\text{d}{v}_x\text{d}{v}_t,(22)$

where r is a zero-mean, unit-variance, delta-correlated, complex Gaussian random function [31], and p and H are

$\begin{array}{c}p\left(v_x,v_t\right)=\sqrt{\frac{\alpha }{\pi }}\exp \left(-\alpha v_x^2\right)\sqrt{\frac{\beta }{\pi }}\exp \left(-\beta v_t^2\right)\\ H\left(\mathit{\rho },t;v_x,v_t\right)=A\exp \left(-\frac{y^2}{4W_y^2}\right)\exp \left(-{\sigma }_x{x}^2\right)\exp \left(-{\sigma }_t{t}^2\right)\\ \exp \left[\text{j}\left(x-\text{j}\alpha \mu t\right)v_x\right]\exp \left[\text{j}\left(t+\text{j}\beta \mu x\right)v_t\right].\end{array}(23)$

The α, β, σ_x, and σ_t relate to the physical twisted space-time beam parameters in Eq. 13 via the relations [19, 35].

$\begin{array}{cc}\hfill \frac{1}{4W_x^2}={\sigma }_x-\frac{\beta {\mu }^2}{2},\hfill & \hfill \frac{1}{4W_t^2}={\sigma }_t-\frac{\alpha {\mu }^2}{2},\hfill \\ \hfill \frac{1}{2{\delta }_x^2}=\frac{\beta {\mu }^2}{4}+\frac{1}{4\alpha },\hfill & \hfill \frac{1}{2{\delta }_t^2}=\frac{\alpha {\mu }^2}{4}+\frac{1}{4\beta }.\hfill \end{array}(24)$

In the simulations, we produced twisted space-time beams with the following parameter values λ_c = 1 μm, W_x = W_y = 2 cm, δ_x = 0.9W_x, W_t = 1 ps, δ_t = 0.9W_t, and $\mu =1/\left({\delta }_x{\delta }_t\right)$—the same as in Figure 1. These parameter values corresponded to α = 3.24 cm², β = 0.81 ps², σ_x = 0.2168 cm⁻², and σ_t = 0.8673 ps⁻². We evaluated Eq. 22 as a matrix-vector product, where the v_x and v_t dimensions were discretized using 64 grid points each, with spacings equal to Δv_x = 0.0645 cm⁻¹ and Δv_t = 0.1291 ps⁻¹, respectively.

3.1.3 Atmospheric turbulence

The index of refraction structure constant and outer scale for the atmospheric turbulence was $C_n^2=10^{-14}{\text{\hspace{0.17em}m}}^{-2/3}$ and L₀ = 10 m, corresponding to the dashed traces in Figure 1. We simulated propagation through this turbulence using the split-step algorithm described in Refs. [70, 77–80]. For N_Fx = 10, 5, 2.5, and 1, we discretized the continuous propagation paths using 4, 5, 9, and 20 equally spaced, statistically independent phase screens generated using the Fourier transform (also known as the spectral) method and augmented with subharmonics [70, 78, 81, 82]. The strength of each phase screen ($C_n^2$, Fried’s parameter r₀, or coherence width ρ₀) was selected such that the discrete-path spherical wave r₀ and scintillation index ${\sigma }_I^2$ matched those of the desired, continuous turbulent path. To capture the change in phase due to turbulence over the light source’s bandwidth, we divided each phase screen by k_c to convert from radians to meters of optical path length (OPL).

Note that we did not simulate the other turbulence conditions reported in Figure 1 due to computational constraints. Accurately simulating turbulence with a given outer scale requires phase screens that have physical dimensions on the order of L₀. Simulating the L₀ = 50 m and 100 m atm would have required grids that were (approximately) 25 and 100 times larger (in numbers of points), respectively, than those used in the L₀ = 10 m simulations (see Section 3.1.1).

3.1.4 Procedure

On each Monte Carlo trial,

1. We generated a twisted space-time beam realization and an instance of atmospheric turbulence as described above.

2. We then Fourier transformed the twisted space-time beam realization to the ω domain using a fast Fourier transform (FFT) computed along the third dimension of U.

3. We propagated U to each of the 4, 5, 9, or 20 (depending on N_Fx) planes using the convolution form of the Fresnel diffraction integral (also known as the angular spectrum propagation method [78, 80]), which we evaluated using FFTs computed along U’s spatial dimensions.

4. In each plane, we converted the atmospheric phase screen in meters of OPL to radians using the ω values along the third dimension of U. We then applied the phase screen to the field and propagated U to the next plane.

5. Upon reaching the observation plane, we Fourier transformed the field back to the t domain using an FFT computed along U’s third dimension.

6. Lastly, we computed the trial irradiance $I\left(x,0,z,t\right)={\left|U\left(x,0,z,t\right)\right|}^2$.

We repeated this procedure 1,000 times.

3.2 Results

Figures 2–4 show the results of the twisted space-time beam simulations. Figures 2, 3—which report the ensemble-averaged irradiances $⟨I\left(x,0,z,t\right)⟩$ after propagating through free space (included as a reference) and atmospheric turbulence, respectively—are organized in the same manner: The top row shows the theoretical $⟨I\left(x,0,z,t\right)⟩$ given in Eq. 19 for Fresnel numbers N_Fx = ∞, 10, 5, and 1, respectively. The bottom (second) row displays the same for the simulated $⟨I\left(x,0,z,t\right)⟩$. The images in Figure 3 are encoded using the same color scales as the corresponding subfigures in Figure 2. Row and column headings have been added to both figures to aid the reader. Lastly, Figure 4 reports the theoretical and simulated ${\widehat{W}}_x/W_x$, ${\widehat{W}}_t/W_t$, and $\widehat{\mu }/\mu$ versus Fresnel number N_Fx. The solid and dashed curves in the figure are the same as those shown in Figure 1B; however, here, we have added the simulated results denoted by the markers ◦ and ⊳. We obtained these results by fitting Gaussian functions to the simulated $⟨I\left(x,0,z,t\right)⟩$.

FIGURE 2

FIGURE 2. Ensemble-averaged irradiance $⟨I\left(x,0,z,t\right)⟩$ free-space results: (A) theory N_Fx = ∞, (B) theory N_Fx = 10, (C) theory N_Fx = 5, (D) theory N_Fx = 1, (E) simulation N_Fx = ∞, (F) simulation N_Fx = 10, (G) simulation N_Fx = 5, and (H) simulation N_Fx = 1.

FIGURE 3

FIGURE 3. Ensemble-averaged irradiance $⟨I\left(x,0,z,t\right)⟩$ turbulence results: (A) theory N_Fx = ∞, (B) theory N_Fx = 10, (C) theory N_Fx = 5, (D) theory N_Fx = 1, (E) simulation N_Fx = ∞, (F) simulation N_Fx = 10, (G) simulation N_Fx = 5, and (H) simulation N_Fx = 1.

FIGURE 4

FIGURE 4. Theory [Eq. 19] and simulation ${\widehat{W}}_x/W_x$, ${\widehat{W}}_t/W_t$, and $\widehat{\mu }/\mu$ versus Fresnel number N_Fx. The symbols ◦ and ⊳ are the results of the simulation.

Inspection of Figure 3 reveals good agreement between simulation and theory in weak to moderately strong atmospheric turbulence [Figures 3B, C, F, and G]. In contrast, the agreement is rather poor in strong turbulence [Figures 3D, H]. This discrepancy is likely caused by the quadratic approximation we used to derive the MCF in Eq. 17 and subsequently $⟨I\left(x,0,z,t\right)⟩$ in Eq. 19. The validity of the quadratic approximation (and the extended Huygens–Fresnel principle, more generally) is suspect in strong turbulence [61–63, 83]. Thus, the disagreement in Figures 3D, H is somewhat expected. The results in Figure 4 are consistent with those in Figure 3—we observe good agreement in weak-to-moderate turbulence and poor agreement in strong turbulence. Although the theoretical relations for ${\widehat{W}}_x/W_x$, ${\widehat{W}}_t/W_t$, and $\widehat{\mu }/\mu$ generally underestimate the effects of turbulence on those parameters, they do accurately predict the trends versus Fresnel number and turbulence strength.

3.3 Experimental verification

Before concluding, we briefly discuss the process for experimentally verifying the theoretical and simulated results presented above. Twisted space-time beam field realizations can be physically synthesized using an apparatus known as a Fourier transform pulse shaper (FTPS) [1, 4, 9, 84–87]. An FTPS consists of two identical gratings separated by a 4f cylindrical lens (CL) system. At the center of the 4f system is a spatial light modulator (SLM). Assuming a pulsed laser beam is input into the FTPS, the first grating-CL-2f system spreads and maps the input beam’s spectrum into physical space at the SLM plane. The SLM modifies the field in the space-frequency (x-ω) domain, which is then transformed back to the space-time domain by the second grating-CL-2f system. Partial coherence manifests by incoherently summing many independent twisted space-time beam realizations.

Turbulence (besides outdoor experiments which are generally uncontrolled) can be controllably generated in the laboratory using several different methods [88]. Of these, phase plate/wheel [89–92] or hot-air [93, 94] techniques are the most germane, and systems employing those methods are easily capable of reproducing the turbulence conditions simulated above.

Lastly, to observe the beam’s behavior in x-t domain, we follow the procedure described in Refs. [1, 5]: The light at the output of the turbulence generator transits a grating-CL-2f system and then is measured by a detector. The detector measures the light’s spatially resolved spectrum averaged over many independent field and turbulence realizations, i.e.,

$S\left(x,z,\omega \right)=⟨{\left|U\left(x,z,\omega \right)\right|}^2⟩.(25)$

Note that this quantity is also referred to as the spectral density [58, 59, 71, 72]. Using Eq. 14, the spectral density relates to the MCF via

$S\left(x,z,\omega \right)=\frac{1}{{\left(2\pi \right)}^2}{\iint }_{-\infty }^{\infty }\mathrm{\Gamma }\left(x,z,t_1,x,z,t_2\right)\exp \left[\text{j}\omega \left(t_1-t_2\right)\right]\text{d}{t}_1\text{d}{t}_2,(26)$

and consequently, the ensemble-averaged irradiance $⟨I\left(x,z,t\right)⟩$ is not directly recoverable. Likely, the easiest course of action is to compare the measured spectral density to its theoretical and simulated counterparts to validate the latter.

4 Conclusion

In this paper, we focused on a recently introduced, partially coherent, space-time-coupled field known as a twisted space-time beam. Twisted space-time beams are similar to traditional twisted Gaussian Schell-model beams; however, instead of being spatially twisted (like the latter), the former possess a stochastic twist which couples their space and time dimensions. Like STOV beams, this spatiotemporal twist imbues twisted space-time beams with transverse (to the direction of propagation) angular momentum.

Generalizing the original research presented in Refs. [19, 20], here, we studied how twisted space-time beams behave as they propagate through atmospheric turbulence. Applying the extended Huygens–Fresnel principle, we derived the MCF for twisted space-time beams after propagating a distance z through atmospheric turbulence of arbitrary strength. From the MCF, we obtained the ensemble-averaged irradiance and quantified the effects of turbulence on beam size, pulse width, and space-time twist. We then simulated twisted space-time beam propagation through atmospheric turbulence to validate our theoretical analysis. The simulated results were found to be in good agreement with theory in weak-to-moderate turbulence. On the other hand, we observed rather poor agreement in strong turbulence, where our theoretical expression for the ensemble-averaged irradiance underestimated the effects of turbulence on the beam size, pulse width, and space-time twist. It did, however, accurately predict the trends of those parameters versus Fresnel number and turbulence strength.

Light with engineered space-time or spatiotemporal coupling is a new and exciting aspect of beam control research, with potential revolutionary uses in optical communications, optical tweezing, and quantum optics. While the free-space propagation characteristics of space-time-coupled beams are generally understood, much less is known about how these beams behave in random media. The results in this paper are a first step toward this goal.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

MH performed the tasks of conceptualization, formal analysis, investigation, methodology, validation, visualization, and writing.

Acknowledgments

The views expressed in this paper are those of the authors and do not reflect the official policy or position of the US Air Force, the Department of Defense, or the US government. MH would like to thank the Air Force Office of Scientific Research (AFOSR) Physical and Biological Sciences Branch (RTB) for supporting this work.

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Kondakci HE, Abouraddy AF. Diffraction-free space-time light sheets. Nat Photon (2017) 11:733–40. doi:10.1038/s41566-017-0028-9

CrossRef Full Text | Google Scholar

2. Yessenov M, Bhaduri B, Kondakci HE, Abouraddy AF. Weaving the rainbow: Space-time optical wave packets. Opt. Photonics News (2019) 30(5):34–41. doi:10.1364/OPN.30.5.000034

CrossRef Full Text | Google Scholar

3. Bhaduri B, Yessenov M, Abouraddy AF. Anomalous refraction of optical spacetime wave packets. Nat Photon (2020) 14:416–21. doi:10.1038/s41566-020-0645-6

CrossRef Full Text | Google Scholar

4. Yessenov M, Hall LA, Schepler KL, Abouraddy AF. Space-time wave packets. Adv Opt Photon (2022) 14:455–570. doi:10.1364/AOP.450016

CrossRef Full Text | Google Scholar

5. Dallaire M, McCarthy N, Piché M. Spatiotemporal bessel beams: Theory and experiments. Opt Express (2009) 17:18148–64. doi:10.1364/OE.17.018148

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Bliokh KY, Nori F. Spatiotemporal vortex beams and angular momentum. Phys Rev A (Coll Park) (2012) 86:033824. doi:10.1103/PhysRevA.86.033824

CrossRef Full Text | Google Scholar

7. Jhajj N, Larkin I, Rosenthal EW, Zahedpour S, Wahlstrand JK, Milchberg HM. Spatiotemporal optical vortices. Phys Rev X (2016) 6:031037. doi:10.1103/PhysRevX.6.031037

CrossRef Full Text | Google Scholar

8. Hancock SW, Zahedpour S, Goffin A, Milchberg HM. Free-space propagation of spatiotemporal optical vortices. Optica (2019) 6:1547–53. doi:10.1364/OPTICA.6.001547

CrossRef Full Text | Google Scholar

9. Chong A, Wan C, Chen J, Zhan Q. Generation of spatiotemporal optical vortices with controllable transverse orbital angular momentum. Nat Photon (2020) 14:350–4. doi:10.1038/s41566-020-0587-z

CrossRef Full Text | Google Scholar

10. Meng X, Hu Y, Wan C, Zhan Q. Optical vortex fields with an arbitrary orbital angular momentum orientation. Opt Lett (2022) 47:4568–71. doi:10.1364/OL.468360

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Wan C, Chen J, Chong A, Zhan Q. Photonic orbital angular momentum with controllable orientation. Natl Sci Rev (2021) 9:nwab149. doi:10.1093/nsr/nwab149

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Yao AM, Padgett MJ. Orbital angular momentum: Origins, behavior and applications. Adv Opt Photon (2011) 3:161–204. doi:10.1364/AOP.3.000161

CrossRef Full Text | Google Scholar

13. Gao D, Ding W, Nieto-Vesperinas M, Ding X, Rahman M, Zhang T, et al. Optical manipulation from the microscale to the nanoscale: Fundamentals, advances and prospects. Light Sci Appl (2017) 6:e17039. doi:10.1038/lsa.2017.39

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Padgett MJ. Orbital angular momentum 25 years on [Invited]. Opt Express (2017) 25:11265–74. doi:10.1364/OE.25.011265

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Shen Y, Wang X, Xie Z, Min C, Fu X, Liu Q, et al. Optical vortices 30 years on: OAM manipulation from topological charge to multiple singularities. Light Sci Appl (2019) 8:90. doi:10.1038/s41377-019-0194-2

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Angelsky OV, Bekshaev AY, Hanson SG, Zenkova CY, Mokhun , Jun Z. Structured light: Ideas and concepts. Front Phys (2020) 8:114. doi:10.3389/fphy.2020.00114

CrossRef Full Text | Google Scholar

17. Yessenov M, Bhaduri B, Kondakci HE, Meem M, Menon R, Abouraddy AF. Non-diffracting broadband incoherent space-time fields. Optica (2019) 6:598–607. doi:10.1364/OPTICA.6.000598

CrossRef Full Text | Google Scholar

18. Mirando A, Zang Y, Zhan Q, Chong A. Generation of spatiotemporal optical vortices with partial temporal coherence. Opt Express (2021) 29:30426–35. doi:10.1364/OE.431882

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Hyde MW. Twisted space-frequency and space-time partially coherent beams. Sci Rep (2020) 10:12443. doi:10.1038/s41598-020-68705-9

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Hyde IVMW. Twisted spatiotemporal optical vortex random fields. IEEE Photon J (2021) 13:1–16. doi:10.1109/JPHOT.2021.3066898

CrossRef Full Text | Google Scholar

21. Ding C, Horoshko D, Korotkova O, Jing C, Qi X, Pan L. Source coherence-induced control of spatiotemporal coherency vortices. Opt Express (2022) 30:19871. doi:10.1364/OE.458666

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Simon R, Mukunda N. Twisted Gaussian schell-model beams. J Opt Soc Am A (1993) 10:95–109. doi:10.1364/JOSAA.10.000095

CrossRef Full Text | Google Scholar

23. Simon R, Sundar K, Mukunda N. Twisted Gaussian Schell-model beams. I. Symmetry structure and normal-mode spectrum. J Opt Soc Am A (1993) 10:2008. doi:10.1364/JOSAA.10.002008

CrossRef Full Text | Google Scholar

24. Sundar K, Simon R, Mukunda N. Twisted Gaussian Schell-model beams. II. Spectrum analysis and propagation characteristics. J Opt Soc Am A (1993) 10:2017–23. doi:10.1364/JOSAA.10.002017

CrossRef Full Text | Google Scholar

25. Simon R, Mukunda N. Twist phase in Gaussian-beam optics. J Opt Soc Am A (1998) 15:2373–82. doi:10.1364/JOSAA.15.002373

CrossRef Full Text | Google Scholar

26. Ambrosini D, Bagini V, Gori F, Santarsiero M. Twisted Gaussian schell-model beams: A superposition model. J Mod Opt (1994) 41:1391–9. doi:10.1080/09500349414551331

CrossRef Full Text | Google Scholar

27. Gori F, Santarsiero M, Borghi R, Vicalvi S. Partially coherent sources with helicoidal modes. J Mod Opt (1998) 45:539–54. doi:10.1080/09500349808231913

CrossRef Full Text | Google Scholar

28. Friberg AT, Tervonen E, Turunen J. Interpretation and experimental demonstration of twisted Gaussian Schell-model beams. J Opt Soc Am A (1994) 11:1818–26. doi:10.1364/JOSAA.11.001818

CrossRef Full Text | Google Scholar

29. Wang H, Peng X, Liu L, Wang F, Cai Y, Ponomarenko SA. Generating bona fide twisted Gaussian Schell-model beams. Opt Lett (2019) 44:3709–12. doi:10.1364/OL.44.003709

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Tian C, Zhu S, Huang H, Cai Y, Li Z. Customizing twisted Schell-model beams. Opt Lett (2020) 45:5880–3. doi:10.1364/OL.405149

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Hyde MW. Stochastic complex transmittance screens for synthesizing general partially coherent sources. J Opt Soc Am A (2020) 37:257–64. doi:10.1364/JOSAA.381772

CrossRef Full Text | Google Scholar

32. Wang H, Peng X, Zhang H, Liu L, Chen Y, Wang F, et al. Experimental synthesis of partially coherent beam with controllable twist phase and measuring its orbital angular momentum. Nanophotonics (2021) 11:689–96. doi:10.1515/nanoph-2021-0432

CrossRef Full Text | Google Scholar

33. Zhang Y, Zhang X, Wang H, Ye Y, Liu L, Chen Y, et al. Generating a twisted Gaussian Schell-model beam with a coherent-mode superposition. Opt Express (2021) 29:41964–74. doi:10.1364/OE.446160

CrossRef Full Text | Google Scholar

34. Gori F, Santarsiero M. Twisted Gaussian Schell-model beams as series of partially coherent modified Bessel-Gauss beams. Opt Lett (2015) 40:1587–90. doi:10.1364/OL.40.001587

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Mei Z, Korotkova O. Random sources for rotating spectral densities. Opt Lett (2017) 42:255–8. doi:10.1364/OL.42.000255

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Borghi R. Twisting partially coherent light. Opt Lett (2018) 43:1627–30. doi:10.1364/OL.43.001627

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Gori F, Santarsiero M. Devising genuine twisted cross-spectral densities. Opt Lett (2018) 43:595–8. doi:10.1364/OL.43.000595

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Serna J, Movilla JM. Orbital angular momentum of partially coherent beams. Opt Lett (2001) 26:405–7. doi:10.1364/OL.26.000405

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Kim SM, Gbur G. Angular momentum conservation in partially coherent wave fields. Phys Rev A (Coll Park) (2012) 86:043814. doi:10.1103/PhysRevA.86.043814

CrossRef Full Text | Google Scholar

40. Gbur GJ. Singular optics. Boca Raton, FL: CRC Press (2016).

Google Scholar

41. Stahl CSD, Gbur G. Twisted vortex Gaussian Schell-model beams. J Opt Soc Am A (2018) 35:1899–906. doi:10.1364/JOSAA.35.001899

CrossRef Full Text | Google Scholar

42. Lin Q, Cai Y. Tensor ABCD law for partially coherent twisted anisotropic Gaussian-Schell model beams. Opt Lett (2002) 27:216–8. doi:10.1364/OL.27.000216

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Cai Y, He S. Propagation of a partially coherent twisted anisotropic Gaussian Schell-model beam in a turbulent atmosphere. Appl Phys Lett (2006) 89:041117. doi:10.1063/1.2236463

CrossRef Full Text | Google Scholar

44. Wang F, Cai Y. Second-order statistics of a twisted Gaussian Schell-model beam in turbulent atmosphere. Opt Express (2010) 18:24661–72. doi:10.1364/OE.18.024661

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wang F, Cai Y, Eyyuboğlu HT, Baykal Y. Twist phase-induced reduction in scintillation of a partially coherent beam in turbulent atmosphere. Opt Lett (2012) 37:184–6. doi:10.1364/OL.37.000184

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Cai Y, Zhu S. Orbital angular moment of a partially coherent beam propagating through an astigmatic ABCD optical system with loss or gain. Opt Lett (2014) 39:1968. doi:10.1364/ol.39.001968

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Charnotskii M. Transverse linear and orbital angular momenta of beam waves and propagation in random media. J Opt (2018) 20:025602. doi:10.1088/2040-8986/aa9f50

CrossRef Full Text | Google Scholar

48. Wang J, Wang H, Zhu S, Li Z. Second-order moments of a twisted Gaussian Schell-model beam in anisotropic turbulence. J Opt Soc Am A (2018) 35:1173–9. doi:10.1364/JOSAA.35.001173

CrossRef Full Text | Google Scholar

49. Zhou M, Fan W, Wu G. Evolution properties of the orbital angular momentum spectrum of twisted Gaussian Schell-model beams in turbulent atmosphere. J Opt Soc Am A (2020) 37:142–8. doi:10.1364/JOSAA.37.000142

CrossRef Full Text | Google Scholar

50. Liu Z, Wan L, Zhou Y, Zhang Y, Zhao D. Progress on studies of beams carrying twist. Photonics (2021) 8:92. doi:10.3390/photonics8040092

CrossRef Full Text | Google Scholar

51. Ponomarenko SA. Classical entanglement of twisted random light propagating through atmospheric turbulence [Invited]. J Opt Soc Am A (2022) 39:C1. C1–C5. doi:10.1364/JOSAA.465410

CrossRef Full Text | Google Scholar

52. Bliokh KY. Spatiotemporal vortex pulses: Angular momenta and spin-orbit interaction. Phys Rev Lett (2021) 126:243601. doi:10.1103/PhysRevLett.126.243601

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Hancock SW, Zahedpour S, Milchberg HM. Mode structure and orbital angular momentum of spatiotemporal optical vortex pulses. Phys Rev Lett (2021) 127:193901. doi:10.1103/PhysRevLett.127.193901

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Lutomirski RF, Yura HT. Propagation of a finite optical beam in an inhomogeneous medium. Appl Opt (1971) 10:1652–8. doi:10.1364/AO.10.001652

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Fante RL, Wolf E. VI wave propagation in random media: A system approach. Opt (Elsevier), Prog Opt (1985) 22:341–98. chap. 6. doi:10.1016/S0079-6638(08)70152-5

CrossRef Full Text | Google Scholar

56. Ishimaru A. Wave propagation and scattering in random media. Piscataway, NJ: IEEE Press (1999). p. 600.

Google Scholar

57. Andrews LC, Phillips RL. Laser beam propagation through random media. 2nd ed. Bellingham, WA: SPIE Press (2005).

Google Scholar

58. Mandel L, Wolf E. Optical coherence and quantum optics. New York, NY: Cambridge University (1995).

Google Scholar

59. Lajunen H, Vahimaa P, Tervo J. Theory of spatially and spectrally partially coherent pulses. J Opt Soc Am A (2005) 22:1536–45. doi:10.1364/josaa.22.001536

CrossRef Full Text | Google Scholar

60. Fante R. Electromagnetic beam propagation in turbulent media. Proc IEEE (1975) 63:1669–92. doi:10.1109/PROC.1975.10035

CrossRef Full Text | Google Scholar

61. Fante R. Electromagnetic beam propagation in turbulent media: An update. Proc IEEE (1980) 68:1424–43. doi:10.1109/PROC.1980.11882

CrossRef Full Text | Google Scholar

62. Fante RL. Two-position, two-frequency mutual-coherence function in turbulence. J Opt Soc Am (1981) 71:1446–51. doi:10.1364/JOSA.71.001446

CrossRef Full Text | Google Scholar

63. Charnotskii M. Extended huygens-fresnel principle and optical waves propagation in turbulence: Discussion. J Opt Soc Am A (2015) 32:1357–65. doi:10.1364/JOSAA.32.001357

CrossRef Full Text | Google Scholar

64.M Abramowitz, and IA Stegun, editors. Handbook of mathematical functions with formulas, graphs, and mathematical tables. Washington, DC: National Bureau of Standards (1964).

Google Scholar

65. Gradshteyn IS, Ryzhik IM. Table of integrals, series, and products. 7th ed. Burlington, MA: Academic Press (2007).

Google Scholar

66. Belafhal A, Chib S, Khannous F, Usman T. Evaluation of integral transforms using special functions with applications to biological tissues. Comp Appl Math (2021) 40:156. doi:10.1007/s40314-021-01542-2

CrossRef Full Text | Google Scholar

67. Young CY, Andrews LC, Ishimaru A. Time-of-arrival fluctuations of a space–time Gaussian pulse in weak optical turbulence: An analytic solution. Appl Opt (1998) 37:7655–60. doi:10.1364/AO.37.007655

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Young CY. Broadening of ultrashort optical pulses in moderate to strong turbulence. In: JC Ricklin, and DG Voelz, editors. Free-space laser communication and laser imaging II, 4821. Seattle, Washington, USA: International Society for Optics and Photonics SPIE (2002). p. 74–81. doi:10.1117/12.452061

CrossRef Full Text | Google Scholar

69. Gbur G, Wolf E. Spreading of partially coherent beams in random media. J Opt Soc Am A (2002) 19:1592–8. doi:10.1364/JOSAA.19.001592

CrossRef Full Text | Google Scholar

70. Gbur G. Partially coherent beam propagation in atmospheric turbulence [Invited]. J Opt Soc Am A (2014) 31:2038–45. doi:10.1364/JOSAA.31.002038

CrossRef Full Text | Google Scholar

71. Korotkova O. Random light beams: Theory and applications. Boca Raton, FL: CRC (2014).

Google Scholar

72. Korotkova O. Theoretical statistical optics. Hackensack, NJ: World Scientific (2021).

Google Scholar

73. Wang F, Liu X, Cai Y. Propagation of partially coherent beam in turbulent atmosphere: A review (invited review). Prog Electromagn Res (2015) 150:123–43. doi:10.2528/PIER15010802

CrossRef Full Text | Google Scholar

74. Wan L, Zhao D. Controllable rotating Gaussian Schell-model beams. Opt Lett (2019) 44:735–8. doi:10.1364/OL.44.000735

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Gori F, Santarsiero M. Devising genuine spatial correlation functions. Opt Lett (2007) 32:3531–3. doi:10.1364/OL.32.003531

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Martínez-Herrero R, Mejías PM, Gori F. Genuine cross-spectral densities and pseudo-modal expansions. Opt Lett (2009) 34:1399–401. doi:10.1364/OL.34.001399

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Martin JM, Flatté SM. Intensity images and statistics from numerical simulation of wave propagation in 3-D random media. Appl Opt (1988) 27:2111–26. doi:10.1364/AO.27.002111

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Schmidt JD. Numerical simulation of optical wave propagation with examples in MATLAB. Bellingham, WA: SPIE Press (2010).

Google Scholar

79. Jarem JM, Banerjee PP. Computational methods for electromagnetic and optical systems. 2nd ed. Boca Raton, FL: CRC Press (2011).

Google Scholar

80. Voelz DG. Computational fourier optics: A matlab tutorial. Bellingham, WA: SPIE (2011).

Google Scholar

81. Lane RG, Glindemann A, Dainty JC. Simulation of a Kolmogorov phase screen. Waves in Random Media (1992) 2:209–24. doi:10.1088/0959-7174/2/3/003

CrossRef Full Text | Google Scholar

82. Frehlich R. Simulation of laser propagation in a turbulent atmosphere. Appl Opt (2000) 39:393–7. doi:10.1364/AO.39.000393

PubMed Abstract | CrossRef Full Text | Google Scholar

83. Fante RL. Some physical insights into beam propagation in strong turbulence. Radio Sci (1980) 15:757–62. doi:10.1029/RS015i004p00757

CrossRef Full Text | Google Scholar

84. Weiner AM. Femtosecond pulse shaping using spatial light modulators. Rev Sci Instrum (2000) 71:1929–60. doi:10.1063/1.1150614

CrossRef Full Text | Google Scholar

85. Ding C, Koivurova M, Turunen J, Setälä T, Friberg AT. Coherence control of pulse trains by spectral phase modulation. J Opt (2017) 19:095501. doi:10.1088/2040-8986/aa7b5e

CrossRef Full Text | Google Scholar

86. Goodman JW. Introduction to fourier optics. 3rd ed. Englewood, CO: Roberts & Company (2005).

Google Scholar

87. Torres-Company V, Lancis J, Andrés P. Space-time analogies in optics. In: E Wolf, editor. Prog. Opt., 56. Elsevier (2011). p. 1–80. chap. 1. doi:10.1016/B978-0-444-53886-4.00001-0

CrossRef Full Text | Google Scholar

88. Jolissaint L. Optical turbulence generators for testing astronomical adaptive optics systems: A review and designer guide. Publ Astron Soc Pac (2006) 118:1205–24. doi:10.1086/507849

CrossRef Full Text | Google Scholar

89. Rodenburg B, Mirhosseini M, Malik M, Magaña-Loaiza OS, Yanakas M, Maher L, et al. Simulating thick atmospheric turbulence in the lab with application to orbital angular momentum communication. New J Phys (2014) 16:033020. doi:10.1088/1367-2630/16/3/033020

CrossRef Full Text | Google Scholar

90. Rickenstorff C, Rodrigo JA, Alieva T. Programmable simulator for beam propagation in turbulent atmosphere. Opt Express (2016) 24:10000–12. doi:10.1364/OE.24.010000

PubMed Abstract | CrossRef Full Text | Google Scholar

91. Dayton D, Spencer M, Hassall A, Rhoadarmer T. Distributed-volume optical disturbance generation in a scaled-laboratory environment using nematic liquid-crystal phase modulators. In: JJ Dolne, and PJ Bones, editors. Unconventional and indirect imaging, image reconstruction, and wavefront sensing 2018, 10772. San Diego, California, USA: International Society for Optics and Photonics SPIE (2018). p. 107720H. doi:10.1117/12.2319968

CrossRef Full Text | Google Scholar

92. Joo JY, Han SG, Lee JH, Rhee HG, Huh J, Lee K, et al. Development and characterization of an atmospheric turbulence simulator using two rotating phase plates. Curr Opt Photon (2022) 6:445–52.

Google Scholar

93. Keskin O, Jolissaint L, Bradley C. Hot-air optical turbulence generator for the testing of adaptive optics systems: Principles and characterization. Appl Opt (2006) 45:4888–97. doi:10.1364/AO.45.004888

PubMed Abstract | CrossRef Full Text | Google Scholar

94. DiComo G, Helle M, Peñano J, Ting A, Schmitt-Sody A, Elle J. Implementation of a long range, distributed-volume, continuously variable turbulence generator. Appl Opt (2016) 55:5192–7. doi:10.1364/AO.55.005192

PubMed Abstract | CrossRef Full Text | Google Scholar