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

Hyde IV, Milo W.

doi:10.3389/fphy.2022.1055401

ORIGINAL RESEARCH article

Front. Phys., 09 January 2023
Sec. Optics and Photonics
Volume 10 - 2022 | https://doi.org/10.3389/fphy.2022.1055401

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

Milo W. Hyde IV*

Department of Engineering Physics, Air Force Institute of Technology, Dayton, OH, United States

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 (ρ, z, ω) = \frac{k}{j 2 π z} \exp (j k z) \iint_{- \infty}^{\infty} U (ρ^{'}, 0, ω) \exp (\frac{j k}{2 z} {|ρ - ρ^{'}|}^{2}) \exp [Ψ (ρ^{'}, 0; ρ, z; ω)] d^{2} ρ^{'}, (1)

where $j = \sqrt{- 1}$ , ω is the radian frequency, k = ω/c is the wavenumber, c is the speed of light, $ρ^{'} = \hat{x} x^{'} + \hat{y} y^{'}$ is the source vector, and $ρ = \hat{x} x + \hat{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 $(ρ^{'}, 0)$ to $(ρ, z)$ 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{gathered} W (ρ_{1}, z, ω_{1}, ρ_{2}, z, ω_{2}) = \frac{k_{1} k_{2}}{{(2 π)}^{2} z^{2}} \exp [j (k_{1} - k_{2}) z] {\iint \iint}_{- \infty}^{\infty} W (ρ_{1}^{'}, 0, ω_{1}, ρ_{2}^{'}, 0, ω_{2}) \\ \exp (\frac{j k_{1}}{2 z} {|ρ_{1} - ρ_{1}^{'}|}^{2}) \exp (- \frac{j k_{2}}{2 z} {|ρ_{2} - ρ_{2}^{'}|}^{2}) \\ ⟨\exp [Ψ (ρ_{1}^{'}, 0; ρ_{1}, z; ω_{1}) + Ψ^{*} (ρ_{2}^{'}, 0; ρ_{2}, z; ω_{2})]⟩ d^{2} ρ_{1}^{'} d^{2} ρ_{2}^{'}, \end{gathered} (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{gathered} ⟨\exp [Ψ (ρ_{1}^{'}, 0; ρ_{1}, z; ω_{1}) + Ψ^{*} (ρ_{2}^{'}, 0; ρ_{2}, z; ω_{2})]⟩ = \exp [- \frac{1}{2} D (ρ_{1}^{'} - ρ_{2}^{'}, 0; ρ_{1} - ρ_{2}, z; ω_{1}, ω_{2})] \\ \approx \exp \{- 2 π^{2} \int_{0}^{z} \int_{0}^{\infty} κ Φ_{n} (κ, ζ) [k_{1}^{2} + k_{2}^{2} - 2 k_{1} k_{2} \exp (- j β κ^{2}) J_{0} (κ R)] d κ d ζ\}, \end{gathered} (3)

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

\begin{aligned} β & = \frac{ζ (z - ζ)}{2 z} (\frac{1}{k_{1}} - \frac{1}{k_{2}}) \\ R & = |\frac{ζ}{z} (ρ_{1} - ρ_{2}) + (1 - \frac{ζ}{z}) (ρ_{1}^{'} - ρ_{2}^{'})| . \end{aligned} (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.,

Γ (ρ_{1}, z, t_{1}, ρ_{2}, z, t_{2}) = \iint_{- \infty}^{\infty} W (ρ_{1}, z, ω_{1}, ρ_{2}, z, ω_{2}) \exp (- j ω_{1} t_{1}) \exp (j ω_{2} t_{2}) d ω_{1} d ω_{2} . (5)

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

\begin{gathered} Γ (ρ_{1}, z, t_{1}, ρ_{2}, z, t_{2}) \\ = \frac{1}{{(2 π)}^{2} c^{2} z^{2}} {\iint \iint}_{- \infty}^{\infty} \iint_{- \infty}^{\infty} ω_{1} ω_{2} W (ρ_{1}^{'}, 0, ω_{1}, ρ_{2}^{'}, 0, ω_{2}) \\ \exp [- \frac{1}{2} D (ρ_{1}^{'} - ρ_{2}^{'}, 0; ρ_{1} - ρ_{2}, z; ω_{1}, ω_{2})] \\ \exp [- j ω_{1} (t_{1} - \frac{z}{c} - \frac{{|ρ_{1} - ρ_{1}^{'}|}^{2}}{2 z c})] \\ \exp [j ω_{2} (t_{2} - \frac{z}{c} - \frac{{|ρ_{2} - ρ_{2}^{'}|}^{2}}{2 z c})] d ω_{1} d ω_{2} d^{2} ρ_{1}^{'} d^{2} ρ_{2}^{'} . \end{gathered} (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{gathered} Γ (ρ_{1}, z, t_{1}, ρ_{2}, z, t_{2}) \approx \frac{1}{λ_{c}^{2} z^{2}} {\iint \iint}_{- \infty}^{\infty} \exp [- j ω_{c} (t_{1} - t_{2} - \frac{{|ρ_{1} - ρ_{1}^{'}|}^{2}}{2 z c} + \frac{{|ρ_{2} - ρ_{2}^{'}|}^{2}}{2 z c})] \\ \iint_{- \infty}^{\infty} W (ρ_{1}^{'}, 0, {\bar{ω}}_{1}, ρ_{2}^{'}, 0, {\bar{ω}}_{2}) \exp [- \frac{1}{2} D (ρ_{1}^{'} - ρ_{2}^{'}, 0; ρ_{1} - ρ_{2}, z; {\bar{ω}}_{1} + ω_{c}, {\bar{ω}}_{2} + ω_{c})] \exp [- j {\bar{ω}}_{1} (t_{1} - \frac{z}{c})] \exp [j {\bar{ω}}_{2} (t_{2} - \frac{z}{c})] d {\bar{ω}}_{1} d {\bar{ω}}_{2} d^{2} ρ_{1}^{'} d^{2} ρ_{2}^{'}, \end{gathered} (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,

Φ_{n} (κ, ζ) = 0.033 C_{n}^{2} (ζ) \frac{\exp (- κ^{2} / κ_{m}^{2})}{{(κ + κ_{0}^{2})}^{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{gathered} D (ρ_{1}^{'} - ρ_{2}^{'}, 0; ρ_{1} - ρ_{2}, z; ω_{1}, ω_{2}) = 0.033 (2 π^{2}) κ_{0}^{- 5 / 3} U (1; \frac{1}{6}; \frac{κ_{0}^{2}}{κ_{m}^{2}}) \\ \times (k_{1}^{2} + k_{2}^{2}) \int_{0}^{z} C_{n}^{2} (ζ) d ζ \\ - 0.033 (4 π^{2}) κ_{0}^{- 5 / 3} k_{1} k_{2} \int_{0}^{z} C_{n}^{2} (ζ) \\ \times \sum_{n = 0}^{\infty} \frac{{(- 1)}^{n}}{n!} {(\frac{κ_{0}^{2} R^{2}}{4})}^{n} U (n + 1; n + \frac{1}{6}; \frac{κ_{0}^{2}}{κ_{m}^{2}}) d ζ, \end{gathered} (9)

where $U (a; c; z)$ 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 (a; c; z)$ . The result, after much analysis, is

\begin{gathered} D (ρ_{1}^{'} - ρ_{2}^{'}, 0; ρ_{1} - ρ_{2}, z; ω_{1}, ω_{2}) \approx 0.7817 κ_{0}^{- 5 / 3} {(k_{1} - k_{2})}^{2} \int_{0}^{z} C_{n}^{2} (ζ) d ζ \\ + 8.7021 κ_{m}^{- 5 / 3} k_{1} k_{2} \int_{0}^{z} C_{n}^{2} (ζ) [{}_{1}F_{1} (- \frac{5}{6}; 1; - \frac{κ_{m}^{2} R^{2}}{4}) - 1] d ζ - 2.3450 κ_{0}^{1 / 3} k_{1} k_{2} \int_{0}^{z} C_{n}^{2} (ζ) R^{2} d ζ . \end{gathered} (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{gathered} D (ρ_{1}^{'} - ρ_{2}^{'}, 0; ρ_{1} - ρ_{2}, z; ω_{1}, ω_{2}) \approx 0.7817 κ_{0}^{- 5 / 3} {(k_{1} - k_{2})}^{2} \int_{0}^{z} C_{n}^{2} (ζ) d ζ \\ + 2.9139 k_{1} k_{2} \int_{0}^{z} C_{n}^{2} (ζ) R^{5 / 3} d ζ - 2.3450 κ_{0}^{1 / 3} k_{1} k_{2} \int_{0}^{z} C_{n}^{2} (ζ) R^{2} d ζ . \end{gathered} (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 $ω_{1} = {\bar{ω}}_{1} + ω_{c}$ and $ω_{2} = {\bar{ω}}_{2} + ω_{c}$ as stipulated in Eq. 7 and noting that $k_{1} k_{2} = ({\bar{k}}_{1} + k_{c}) ({\bar{k}}_{2} + k_{c}) \approx k_{c}^{2}$ , we arrive at the final result

\begin{gathered} D (ρ_{1}^{'} - ρ_{2}^{'}, 0; ρ_{1} - ρ_{2}, z; {\bar{ω}}_{1} + ω_{c}, {\bar{ω}}_{2} + ω_{c}) \\ \approx \frac{0.7817 C_{n}^{2} z κ_{0}^{- 5 / 3}}{c^{2}} {({\bar{ω}}_{1} - {\bar{ω}}_{2})}^{2} + 1.0930 C_{n}^{2} z (1 - 0.7152 κ_{0}^{1 / 3}) k_{c}^{2} \\ [{|ρ_{1} - ρ_{2}|}^{2} + {|ρ_{1}^{'} - ρ_{2}^{'}|}^{2} + (ρ_{1} - ρ_{2}) \cdot (ρ_{1}^{'} - ρ_{2}^{'})] \\ = 2 \frac{a_{ω}}{c^{2}} {({\bar{ω}}_{1} - {\bar{ω}}_{2})}^{2} + 2 a_{s} k_{c}^{2} [{|ρ_{1} - ρ_{2}|}^{2} + {|ρ_{1}^{'} - ρ_{2}^{'}|}^{2} + (ρ_{1} - ρ_{2}) \cdot (ρ_{1}^{'} - ρ_{2}^{'})] . \end{gathered} (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{gathered} Γ (ρ_{1}, t_{1}, ρ_{2}, t_{2}) = A^{2} \exp (- \frac{y_{1}^{2} + y_{2}^{2}}{4 W_{y}^{2}}) \exp (- \frac{x_{1}^{2} + x_{2}^{2}}{4 W_{x}^{2}}) \exp [- \frac{{(x_{1} - x_{2})}^{2}}{2 δ_{x}^{2}}] \\ \exp (- \frac{t_{1}^{2} + t_{2}^{2}}{4 W_{t}^{2}}) \exp [- \frac{{(t_{1} - t_{2})}^{2}}{2 δ_{t}^{2}}] \exp [j μ (x_{1} t_{2} - x_{2} t_{1})] \exp [- j ω_{c} (t_{1} - t_{2})]; \end{gathered} (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 $|μ| δ_{t} δ_{x} \leq 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 $|μ| δ_{t} δ_{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 (ρ_{1}, ω_{1}, ρ_{2}, ω_{2}) = \frac{1}{{(2 π)}^{2}} \iint_{- \infty}^{\infty} Γ (ρ_{1}, t_{1}, ρ_{2}, t_{2}) \exp (j ω_{1} t_{1}) \exp (- j ω_{2} t_{2}) d t_{1} d t_{2} . (14)

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

\begin{gathered} W (ρ_{1}, ω_{1}, ρ_{2}, ω_{2}) = \frac{A^{2}}{4 π Ω_{t}} \exp (- \frac{y_{1}^{2} + y_{2}^{2}}{4 W_{y}^{2}}) \exp (- \frac{x_{1}^{2} + x_{2}^{2}}{4 {\tilde{W}}_{x}^{2}}) \exp [- \frac{{(x_{1} - x_{2})}^{2}}{2 {\tilde{δ}}_{x}^{2}}] \\ \exp (- \frac{{\bar{ω}}_{1}^{2} + {\bar{ω}}_{2}^{2}}{4 W_{ω}^{2}}) \exp [- \frac{{({\bar{ω}}_{1} - {\bar{ω}}_{2})}^{2}}{2 δ_{ω}^{2}}] \exp \{\frac{μ}{δ_{ω}^{2}} [{\bar{ω}}_{2} (1 + \frac{1}{2 γ_{t}^{2}}) - {\bar{ω}}_{1}] x_{1}\} \\ \exp \{\frac{μ}{δ_{ω}^{2}} [{\bar{ω}}_{1} (1 + \frac{1}{2 γ_{t}^{2}}) - {\bar{ω}}_{2}] x_{2}\}, \end{gathered} (15)

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

\begin{gathered} Ω_{t}^{2} = {(\frac{1}{4 W_{t}^{2}} + \frac{1}{2 δ_{t}^{2}})}^{2} - {(\frac{1}{2 δ_{t}^{2}})}^{2}, \\ \frac{1}{{\tilde{W}}_{x}^{2}} = \frac{1}{W_{x}^{2}} + \frac{μ^{2}}{W_{ω}^{2}}, \\ \frac{1}{{\tilde{δ}}_{x}^{2}} = \frac{1}{δ_{x}^{2}} + \frac{μ^{2}}{δ_{ω}^{2}} . \end{gathered} (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{gathered} Γ (x_{1}, 0, z, t_{1}, x_{2}, 0, z, t_{2}) = A^{2} \frac{W_{t}}{{\tilde{W}}_{t}} \frac{N_{F x} N_{F y}}{\sqrt{Δ_{x} Δ_{y}}} \exp [- j ω_{c} ({\bar{t}}_{1} - {\bar{t}}_{2})] \\ \exp (- \frac{N_{F x}^{2}}{Δ_{x}} \frac{x_{1}^{2} + x_{2}^{2}}{4 W_{x}^{2}}) \exp [- \frac{{\bar{t}}_{1}^{2} + {\bar{t}}_{2}^{2}}{4 {(W_{t}^{eff})}^{2}}] \exp [μ \frac{N_{F x}}{Δ_{x}} \frac{W_{t}^{2}}{{\tilde{W}}_{t}^{2}} (x_{1} {\bar{t}}_{2} + x_{2} {\bar{t}}_{1})] \\ \exp [\frac{j k_{c}}{2 z} (1 - \frac{N_{F x}^{2}}{Δ_{x}} + a_{s} k_{c}^{2} \frac{4 W_{x}^{2}}{Δ_{x}}) (x_{1}^{2} - x_{2}^{2})] \exp [j μ^{2} \frac{N_{F x}^{2}}{Δ_{x}} \frac{W_{t}^{2}}{{\tilde{W}}_{t}^{2}} ({\bar{t}}_{1}^{2} - {\bar{t}}_{2}^{2})] \\ \exp [j μ \frac{N_{F x}^{2}}{Δ_{x}} \frac{W_{t}^{2}}{{\tilde{W}}_{t}^{2}} (x_{1} {\bar{t}}_{2} - x_{2} {\bar{t}}_{1})] \exp [- j μ \frac{N_{F x}^{2}}{Δ_{x}} \frac{a_{ω}}{c^{2} {\tilde{W}}_{t}^{2}} (x_{1} + x_{2}) ({\bar{t}}_{1} - {\bar{t}}_{2})] \\ \exp [- j μ \frac{2 W_{x}^{2}}{Δ_{x}} \frac{W_{t}^{2}}{{\tilde{W}}_{t}^{2}} a_{s} k_{c}^{2} (x_{1} - x_{2}) ({\bar{t}}_{1} + {\bar{t}}_{2})] \\ \exp [- \frac{{(x_{1} - x_{2})}^{2}}{2 {(δ_{x}^{eff})}^{2}}] \exp [- \frac{{({\bar{t}}_{1} - {\bar{t}}_{2})}^{2}}{2 {(δ_{t}^{eff})}^{2}}] \exp [- \frac{(x_{1} - x_{2}) ({\bar{t}}_{1} - {\bar{t}}_{2})}{2 {(δ_{x t}^{eff})}^{2}}] . \end{gathered} (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_{F x, y} = 2 k_{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} + 2 a_{ω} / c^{2}$ ; and

\begin{gathered} Δ_{x} = 1 + 4 \frac{W_{x}^{2}}{δ_{x}^{2}} + N_{F x}^{2} + 8 W_{x}^{2} (a_{s} k_{c}^{2} + μ^{2} \frac{a_{ω}}{c^{2}} \frac{W_{t}^{2}}{{\tilde{W}}_{t}^{2}}); \\ Δ_{y} = 1 + N_{F y}^{2} + 8 W_{y}^{2} a_{s} k_{c}^{2}; \\ \frac{1}{4 {(W_{t}^{eff})}^{2}} = \frac{1}{4 {\tilde{W}}_{t}^{2}} + μ^{2} \frac{W_{x}^{2}}{Δ_{x}} \frac{W_{t}^{4}}{{\tilde{W}}_{t}^{4}}; \\ \frac{1}{2 {(δ_{t}^{eff})}^{2}} = \frac{1}{2 δ_{t}^{2}} + μ^{2} \frac{W_{x}^{2}}{2 Δ_{x}} (Δ_{x} - N_{F x}^{2} - \frac{W_{t}^{4}}{{\tilde{W}}_{t}^{4}}) + \frac{a_{ω}}{c^{2}} \frac{1}{4 W_{t}^{2} {\tilde{W}}_{t}^{2}}; \\ \frac{1}{2 {(δ_{x}^{eff})}^{2}} = \frac{N_{F x}^{2}}{Δ_{x}} \frac{1}{2 δ_{x}^{2}} + a_{s} k_{c}^{2} (1 + 2 \frac{N_{F x}^{2}}{Δ_{x}}) - {(a_{s} k_{c}^{2})}^{2} \frac{2 W_{x}^{2}}{Δ_{x}^{2}} + μ^{2} \frac{a_{ω}}{c^{2}} \frac{N_{F x}^{2}}{Δ_{x}} \frac{W_{t}^{2}}{{\tilde{W}}_{t}^{2}}; \\ \frac{1}{2 {(δ_{x t}^{eff})}^{2}} = μ \frac{N_{F x}}{2 Δ_{x}} (Δ_{x} - N_{F x}^{2} - \frac{W_{t}^{2}}{{\tilde{W}}_{t}^{2}} + a_{s} k_{c}^{2} 4 W_{x}^{2}) . \end{gathered} (18)

Eq. 17 is organized so that the terms can be physically interpreted: Starting at the top and ignoring the carrier $\exp [- j ω_{c} ({\bar{t}}_{1} - {\bar{t}}_{2})]$ , 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{gathered} ⟨I (x, 0, z, t)⟩ = Γ (x, 0, z, t, x, 0, z, t) \\ = A^{2} \frac{W_{t}}{{\tilde{W}}_{t}} \frac{N_{F x} N_{F y}}{\sqrt{Δ_{x} Δ_{y}}} \exp (- \frac{N_{F x}^{2}}{Δ_{x}} \frac{x^{2}}{2 W_{x}^{2}}) \exp [- (1 + μ^{2} \frac{4 W_{x}^{2} W_{t}^{4}}{Δ_{x} {\tilde{W}}_{t}^{2}}) \frac{{\bar{t}}^{2}}{2 {\tilde{W}}_{t}^{2}}] \exp (μ \frac{2 N_{F x}}{Δ_{x}} \frac{W_{t}^{2}}{{\tilde{W}}_{t}^{2}} x \bar{t}) \\ = {\hat{A}}^{2} \exp (- \frac{x^{2}}{2 {\hat{W}}_{x}^{2}}) \exp (- \frac{{\bar{t}}^{2}}{2 {\hat{W}}_{t}^{2}}) \exp (\hat{μ} x \bar{t}) . \end{gathered} (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 ${\hat{W}}_{x}$ , ${\hat{W}}_{t}$ , and $\hat{μ}$ versus Fresnel number and turbulence strength. Figure 1 shows these curves: (A) plots ${\hat{W}}_{x} / W_{x}$ , ${\hat{W}}_{t} / W_{t}$ , and $\hat{μ} / μ$ 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 m^{- 2 / 3}$ ) and atmospheric turbulence ( $C_{n}^{2} = 1 0^{- 14} 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.,

σ_{I}^{2} = 0.5 C_{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 ${\hat{W}}_{x} / W_{x}$ , ${\hat{W}}_{t} / W_{t}$ , and $\hat{μ} / μ$ 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 $μ = 1 / (δ_{x} δ_{t})$ .

FIGURE 1

FIGURE 1. (A) ${\hat{W}}_{x} / W_{x}$ , ${\hat{W}}_{t} / W_{t}$ , and $\hat{μ} / μ$ from Eq. 19 versus Fresnel number N_Fx. (B) Zoomed-in view of ${\hat{W}}_{x} / W_{x}$ , ${\hat{W}}_{t} / W_{t}$ , and $\hat{μ} / μ$ over the boxed region in (A). The solid, dashed, dashed-dotted, and dotted traces show the results in free space $(C_{n}^{2} = 0 m^{- 2 / 3})$ and atmospheric turbulence ( $C_{n}^{2} = 1 0^{- 14} 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 $σ_{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 $μ = 1 / (δ_{x} δ_{t})$ .

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 ${\hat{W}}_{x} \approx W_{x}$ , ${\hat{W}}_{t} \approx W_{t}$ , and $\hat{μ} < μ / 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 ${\hat{W}}_{x}$ , the pulse width also changes in this region because of spatiotemporal coupling. Beginning around N_Fx ≈ 10, ${\hat{W}}_{t}$ starts to contract (shorten) and continues to do so until N_Fx ≈ 1. This shortening of ${\hat{W}}_{t}$ is met by an increase in $\hat{μ}$ . When considered together, the result is a beam that rotates in the x-t (or x-z) plane—the beam effectively “trades” ${\hat{W}}_{t}$ to do so. Lastly, for N_Fx < 1, ${\hat{W}}_{x}$ continues to grow larger due to diffraction, ${\hat{W}}_{t}$ asymptotes (the pulse width stops contracting), and $\hat{μ}$ 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 ${\hat{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 ${\hat{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 ${\hat{W}}_{t}$ , namely,

\frac{1}{{\hat{W}}_{t}} = \frac{1}{{\tilde{W}}_{t}} (1 + μ^{2} \frac{4 W_{x}^{2} W_{t}^{4}}{Δ_{x} {\tilde{W}}_{t}^{2}}) . (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 ${\hat{W}}_{t} \sim {\tilde{W}}_{t} = W_{t} + 2 a_{ω} / 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 $\hat{μ}$ decays much more rapidly in turbulence than in free space. Examining the mathematical relation for $\hat{μ}$ reveals that it approaches zero like $\hat{μ} \sim z^{- 3}$ in turbulence (vice $\hat{μ} \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 (x, 0, z, t)⟩$ 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 (ρ, t) = \iint_{- \infty}^{\infty} r (v_{x}, v_{t}) \sqrt{\frac{1}{2} p (v_{x}, v_{t})} H (ρ, t; v_{x}, v_{t}) d v_{x} 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{gathered} p (v_{x}, v_{t}) = \sqrt{\frac{α}{π}} \exp (- α v_{x}^{2}) \sqrt{\frac{β}{π}} \exp (- β v_{t}^{2}) \\ H (ρ, t; v_{x}, v_{t}) = A \exp (- \frac{y^{2}}{4 W_{y}^{2}}) \exp (- σ_{x} x^{2}) \exp (- σ_{t} t^{2}) \\ \exp [j (x - j α μ t) v_{x}] \exp [j (t + j β μ x) v_{t}] . \end{gathered} (23)

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

\begin{matrix} \frac{1}{4 W_{x}^{2}} = σ_{x} - \frac{β μ^{2}}{2}, & \frac{1}{4 W_{t}^{2}} = σ_{t} - \frac{α μ^{2}}{2}, \\ \frac{1}{2 δ_{x}^{2}} = \frac{β μ^{2}}{4} + \frac{1}{4 α}, & \frac{1}{2 δ_{t}^{2}} = \frac{α μ^{2}}{4} + \frac{1}{4 β} . \end{matrix} (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 $μ = 1 / (δ_{x} δ_{t})$ —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} = 1 0^{- 14} 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 $σ_{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 (x, 0, z, t) = {|U (x, 0, z, t)|}^{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 (x, 0, z, t)⟩$ 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 (x, 0, z, t)⟩$ 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 (x, 0, z, t)⟩$ . 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 ${\hat{W}}_{x} / W_{x}$ , ${\hat{W}}_{t} / W_{t}$ , and $\hat{μ} / μ$ 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 (x, 0, z, t)⟩$ .

FIGURE 2

FIGURE 2. Ensemble-averaged irradiance $⟨I (x, 0, z, t)⟩$ 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 (x, 0, z, t)⟩$ 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 ${\hat{W}}_{x} / W_{x}$ , ${\hat{W}}_{t} / W_{t}$ , and $\hat{μ} / μ$ 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 (x, 0, z, t)⟩$ 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 ${\hat{W}}_{x} / W_{x}$ , ${\hat{W}}_{t} / W_{t}$ , and $\hat{μ} / μ$ 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 (x, z, ω) = ⟨{|U (x, z, ω)|}^{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 (x, z, ω) = \frac{1}{{(2 π)}^{2}} \iint_{- \infty}^{\infty} Γ (x, z, t_{1}, x, z, t_{2}) \exp [j ω (t_{1} - t_{2})] d t_{1} d t_{2}, (26)

and consequently, the ensemble-averaged irradiance $⟨I (x, z, t)⟩$ 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

ORIGINAL RESEARCH article

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

1 Introduction

2 Theory

2.1 Extended Huygens–Fresnel principle

2.2 Approximate two-point, spherical WSF D

2.3 CSD function of a twisted space-time beam

2.4 MCF of twisted space-time beam in atmospheric turbulence

2.5 Average irradiance and physical discussion

3 Validation

3.1 Simulation setup

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

3.1.2 Generating twisted space-time fields

3.1.3 Atmospheric turbulence

3.1.4 Procedure

3.2 Results

3.3 Experimental verification

4 Conclusion

Data availability statement

Author contributions

Acknowledgments

Conflict of interest

Publisher’s note

References

This article is part of the Research Topic

People also looked at