Twisted Spatiotemporal Optical Vortex Random Fields

We present twisted spatiotemporal optical vortex (STOV) beams, which are partially coherent light sources that possess a coherent optical vortex and a random twist coupling their space and time dimensions. These beams have controllable partial coherence and transverse orbital angular momentum (OAM), which distinguishes them from the more common spatial vortex and twisted beams (known to carry longitudinal OAM) in the literature and should ultimately make them useful in applications such as optical communications and optical tweezing. We present the mathematical analysis of twisted STOV beams, deriving the mutual coherence function and linear and angular momentum densities. We simulate the synthesis of twisted STOV beams and investigate their free-space propagation characteristics. We discuss how to physically generate twisted STOV fields and lastly conclude with a summary and brief discussion of future research.

To date, much of the OAM beam research has focused on fields with wavefront twists or vortices coupling their transverse (to the direction of propagation) spatial dimensions. This results in beams with axial or longitudinal (parallel to the propagation direction) OAM. As a result, there has been significant interest in generating optical fields which possess transverse angular momentum. It has been known for decades that evanescent waves excited by total internal reflection and tightly focused Gaussian beams carry transverse spin angular momentum (SAM), also called photonic wheels [30]- [32]. However, only in the last few years have fields with transverse SAM been demonstrated for optical trapping, particle manipulation, optical field probing/reconstruction, and spin-direction locking [18], [30]- [36]. More information on light that possesses transverse SAM and potential applications can be found in Refs. [18], [30]- [32], [37].
Optical beams that carry transverse OAM are a new phenomenon. In series of recent papers, spatiotemporal optical vortex (STOV) and twisted space-time light beams were introduced [38]- [44]. In the former case, an optical vortex couples the space and time dimensions of a coherent light source; whereas in the latter, the space-time coupling (the twist) is random, such that the light beam is partially coherent. The difference in OAM direction between vortex and twisted spatiotemporal beams and their more common spatial (longitudinal OAM) counterparts discussed above creates possibilities for STOV and twisted space-time fields to be applied in optical tweezing, atomic optics, optical communications, and laser structuring in novel ways.
In this article, we combine STOV and twisted space-time beams to produce twisted STOV fields. These partially coherent light beams possess both a random twist and a coherent vortex coupling their space and time dimensions. This combination of controllable partial coherence and transverse OAM should make these beams useful in existing OAM beam applications.
We begin with the mathematical derivation of twisted STOV beams and obtain expressions for their mutual coherence function (MCF) and linear and angular momentum densities. We then simulate the synthesis of twisted STOV beams and examine their free-space propagation properties. We lastly conclude with a brief summary and a discussion of future work.

Source-Plane Twisted STOV MCF
To find the MCF of a twisted STOV beam, we begin with the superposition rule for partially coherent fields originally derived by Gori and Santarsiero [45]: where, in general, p is any real positive function and H is a kernel. As is customary for space-time coupled beams, we ignore the beam's spatial distribution in the y direction. If we remove the vortex from a twisted STOV beam, the MCF should simplify to that of a twisted space-time partially coherent field described in Ref. [44]. Thus, we should expect that the p and H in Eq. (1) are very similar to those of twisted space-time beams. We therefore let p and H equal where σ x , σ t , α, and β are related to the beam's spatial and temporal pulse and correlation widths (W x , W t , δ x , and δ t , respectively), and μ is the twist parameter subject to the constraint |μ|δ x δ t ≤ 1 [46]. With the exception of the complex function τ (described more below), the p and H in Eq. (2) are the same as those in Ref. [44]. We can therefore substitute them into Eq. (1) and evaluate the integrals to yield The p and H parameters (σ x , σ t , α, and β) relate to physical beam parameters W x , W t , δ x , and δ t through the following system of nonlinear equations: which unfortunately, for synthesis of twisted space-time beams, cannot be inverted. The MCF in Eq. (3) is organized such that the first four terms correspond to the beam's deterministic spatial and temporal (in general, complex) shapes, while the latter three define the spatial and temporal correlation properties of the random beam. STOV beams, as discussed in Refs. [39]- [43], are fully coherent, deterministic optical fields and take a form similar to that of Laguerre-Gauss laser modes [39], [40], [42]. Thus, the spatiotemporal vortex should be completely described in the first four terms of the MCF in Eq. (3). In addition, those terms should assume the form of a Laguerre-Gauss mode [42], [47], [48], namely, where ρ and φ are the spatial-temporal radial and azimuthal coordinates, ω c is the mean, or carrier frequency of the light beam, W is the beam's spatial-temporal radius, and is the STOV order. In Eq. (5), we have assumed a Laguerre-Gauss mode with radial index equal to zero. From here, it is a rather simple matter of finding ρ, φ, and τ ; they are where sgn(x ) is the signum function. The twisted STOV MCF becomes where ρ =xxW/(2W x ) +ttW/(2W t ) andx ×t =ŷ, ξ = 4μW x W t /W 2 is the overall space-time twist, is the spatial-temporal correlation radius. In addition to the relations in Eq. (4), ασ x = βσ t .

Propagation of the Twisted STOV MCF
Not surprisingly, Eq. (7) is very similar in form to the source-plane (z = 0) cross-spectral density function for spatially twisted Laguerre-Gaussian Schell-model beams [49], [50]. Despite this similarity, twisted STOV beams propagate differently than traditional spatially twisted partially coherent fields because of spatiotemporal coupling. Assuming the field is relatively narrowband, i.e., ω c max(1/W t , 1/δ t ), the twisted STOV MCF at any transverse plane z > 0 in free space can be found by evaluating where c is the speed of light and λ c is the beam's carrier wavelength. Substituting Eq. (7) into Eq. (8) yields where Following approaches similar to those employed in Refs. [49], [51], it might be possible to evaluate Eq. (9) analytically. The general procedure would be to expand the th-order polynomials using the binomial theorem and then evaluate the resulting powers-times-Gaussian integrals. Ultimately, the MCF would be expressed in the form of multidimensional sums. Computationally, there is little difference evaluating multidimensional sums vice integrals; thus, we choose the latter.
In the following section, we discuss the details of this computation. We also explain how to synthesize realizations of twisted STOV beams and discuss the particulars of the wave-optics simulations. The theoretical MCF and simulation results are then presented in Section IV.

Linear and Angular Momentum Densities
Before proceeding to the computational details, it is worth discussing the OAM of twisted STOV beams. We begin with the linear P and angular L momentum densities for coherent paraxial beams propagating in the z direction originally derived by Allen et al. [1]: where ε 0 is the vacuum permittivity and r =ρρ +ẑz is the position vector. The first term of P, after neglecting the relatively weak ∂U/∂z, is the transverse linear momentum density. For beams with an x-y (spatial) vortex or twist, this term, via the cross-product with r and subsequent integration over the beam's cross section, gives rise to a z-directed total OAM [1], [3], [9], [51]- [53]. For twisted STOV beams, the mathematics is similar; however, as we show, the total OAM is y directed.
Twisted STOV beams are partially coherent, narrowband fields, and therefore, we are interested in the ensemble averages of P and L, such that where ∇ 2 operates on the subscript "2" coordinates of the MCF [9], [52], [53]. Using Eqs. (6)-(9), the del operator in spatial-temporal cylindrical coordinates is where the t and z directions are parallel viat , such that W z = cW t (note, in general, where v g is the group velocity of the pulse) and the twisted STOV MCF is assumed invariant in the y direction. As a result of this last assumption, the position vector r = ρ =xxW/(2W x ) +ẑzW/(2W z ). Substituting Eqs. (7) (the source-plane MCF) and (12) into Eq. (11) and simplifying yields the source-plane linear momentum density: We can now derive the angular momentum density L by taking the cross product of r and P , viz., By noting that xW/(2W x ) = ρ cos φ and zW/(2W z ) = ρ sin φ, L simplifies to Integration over the x and z (or t ) dimensions of the beam yields the total OAM: Like the source-plane MCF given in Eq. (7), this expression for the total OAM is very similar in form to that derived in Refs. [49], [51] for spatially twisted vortex beams. The critical difference being direction-longitudinal or axial OAM in Refs. [49], [51], transverse OAM here. It is well known that beams which carry longitudinal OAM rotate in the transverse plane as they propagate [6], [54], [55]. Similarly, twisted STOV beams-carrying transverse OAM-rotate in a plane parallel to the propagation direction. We discuss this in more detail in Section IV.

Computing the Twisted STOV MCF At Any Plane z > 0
We begin this section with the numerical calculation of the MCF in Eq. (9). The integrand in Eq. (9) is expressed as the product of a six-dimensional function (lines 2 and 3) h(x 1 , t 1 , x 2 , t 2 ; x 1 , x 2 ) and a two-dimensional (2D) function χ (x 1 , x 2 ) (line 4). The subsequent integrals over x 1 and x 2 form a 2D superposition integral. Equation (9), in its entirety, can be evaluated as the Hadamard product of a matrix-vector product, namely, where x 1 and x 2 are the grid spacings in the source x 1 -x 2 plane;t is a vector representing the "t terms" (line 1); H is a matrix corresponding to h with the observation (x 1 , t 1 , x 2 , and t 2 ) and source coordinates along the rows and columns, respectively; and lastly, x is a vector representing χ . The MCF in Eq. (17) is a vector which should be reshaped into a more physically meaningful form.
The main challenge with evaluating Eq. (9) in this manner is the size of H. Indeed, H becomes prohibitive for even small computational grids, and it is therefore impractical to find at every combination of x 1 , t 1 , x 2 , and t 2 . Fortunately, the MCF's behavior as a function of z can be determined from 2D planar cuts through the four-dimensional (4D) MCF. This reduces H to a more manageable size. In the theoretical MCF results presented in the next section, we compute Eq. (9), for several propagation distances, under the following conditions: r Letting x 1 = x 2 and t 1 = t 2 , we find the ensemble-averaged, or pulse-averaged intensity r Letting x 2 = 0 and t 1 = z/c (i.e.,t 1 ≈ 0), we find (x 1 , 0, z/c, t 2 , z), which shows the twisted STOV beam's spatiotemporal coupling. The particular twisted STOV beam parameter values (W x , W t , δ x , δ t , μ, and ), computational grid sizes and spacings, and propagation distances are presented at the end of this section.

Generating Twisted STOV Realizations
We now turn to generating stochastic realizations of twisted STOV beams. The goal is to generate optical fields which are sample functions drawn from the set of all such functions whose ensembleaveraged autocorrelation is given by the MCF in Eq. (7). Starting with a vector of zero-mean, unit-variance, independent, complex Gaussian random numbers r, a twisted STOV field instance can be generated by evaluating where p is given in Eq.
Since r is delta correlated, i.e., , the 2D integral in Eq. (18) can be evaluated as the product of two one-dimensional integrals. In discrete form, Eq. (18) becomes where v x and v t are the spacings in the v x and v t directions, τ is a vector representing the functions on line 1 of Eq. (18), and V x and V t are matrices corresponding to the Fourier-like kernels in Eq. (18). The other symbols are vectors corresponding to r and p. Like , U is a vector which should be reshaped into a matrix to physically represent the optical field. In contrast to Eq. (18), where a linear operation is performed on a matrix of Gaussian random numbers, Eq. (19) contains the product of two independent normal random vectors. Therefore, the U generated using Eq. (19) are not Gaussian distributed. Nevertheless, Eq. (19) produces U realizations with the desired twisted STOV MCF.

Simulation Setup
With Eqs. (17) and (19), we can now examine the propagation behaviors of twisted STOV beams. We considered three circular (i.e., W x = cW t ) twisted STOV beams in this analysis: an "incoherent" case where δ x = 3W x /4 and δ t = 3W t /4, a "partially coherent" case where δ x = W x and δ t = W t , and a  (4), |μ|δ x δ t ≤ 1, and ασ x = βσ t ], the twist μ was forced to vary from case to case. For this study, we would have preferred to hold μ constant; however, we were able to saturate the twists, such that |μ|δ x δ t = 1 [57]. In addition, we set μ < 0 so that the twist OAM component was in the same direction as that of the vortex [see Eqs. (15) and (16)]. All parameter values are reported in Table 1.
We propagated these three beams to five z locations, i.e., z = 70.59 m, 141.18 m, 282.35 m, 564.70 m, and 1129.41 m, corresponding to Fresnel numbers N F = 4πW 2 x /(λ c z) = 10, 5, 2.5, 1.25, and 0.625, to examine how the beams evolved from the near to far zones. We found the theoretical MCFs by computing Eq. (17) as described above.
For the twisted STOV field instances U , generated via Eq. (19), we first transformed U to the x-ω domain by performing a fast Fourier transform (FFT) along the t dimension of U . We then propagated U (x, ω) to the z locations listed above by evaluating the Fresnel integral [58]- [60] using a FFT computed along U 's x dimension. Lastly, we transformed U (x, ω, z) back to the x-t domain using a FFT computed along U 's ω dimension resulting in U (x, t , z). From 15000 statistically independent U (x, t , z), we computed the pulse-averaged intensity I(x, t , z) and (x 1 , 0, z/c, t 2 , z) to compare to Eq. (17).
In both the source and observation planes, we used grids that were N = 512 points per side with spacings equal to x 1 = x 2 = x = min(δ x )/10 and x 1 = x 2 = x = λ c z/(N x ), respectively. For the t dimension, the spacings were t = min(δ t )/10 in both planes. Lastly, in Eq. (19), the spacings in the v x -v t domain were v x = 2π/(N x ) and v t = 2π/(N t ).
We have included the MATLAB.m files required to execute the simulations as supplementary materials to this paper. In the next section, we present and discuss the results.

Twisted STOV Fields
We begin this section with twisted STOV field realizations shown in Fig. 1. These are included to show what twisted STOV fields physically look like. Each row of the figure displays the magnitude and phase of an example source-plane U (x, t ) from the incoherent (I), partially coherent (PC), and coherent (C) beam cases described above. Row and column headings are included for the reader's convenience.
Examining the |U | images, the similarity of these beams, for lack of a better comparison, to Swiss cheese is striking. Of course, the "holes" are caused by the phase vortices evident in the arg(U ) images. In each of these plots, one can clearly discern the second-order STOV at the origin. Being deterministic, this phase dislocation is present in every twisted STOV realization.
The random twist arises from the first-order vortices that surround the on-axis STOV. The locations of these "twist" vortices change from realization to realization, and clearly, the vortex density decreases as the field becomes more coherent. In addition, all the twist vortices in the arg(U ) images have the same topological charge (positive in Fig. 1), and therefore, connect to oppositely (negatively) charged vortices at infinity [9], [61], [62]. This produces the web-like patterns in the arg(U ) plots.
All twist vortices have the same charge because the twist is saturated in the I, PC, and C cases. In the opposite extreme, setting μ = 0 and = 0, the twisted STOV MCF in Eq. (7) simplifies to that of a Gaussian Schell-model (GSM) pulsed beam [63]- [65]. Fig. 2 shows GSM pulsed beam realizations for the I, PC, and C beam cases. These instances were generated in the same manner as the twisted STOV realizations described above. GSM pulsed beam realizations are space-time speckle fields, and like spatial speckle fields, contain roughly equal numbers of positively and negatively charged phase vortices [9], [61], [62]. This means, as |μ|δ x δ t < 1 and ultimately approaches 0, twisted STOV field realizations (excluding the = 0 deterministic on-axis vortex) should evolve into traditional speckle fields, and more and more negatively charged vortices (in our case) should appear in the arg(U ) images.     These results prove that we are indeed producing twisted STOV field realizations with the proper statistics.

Propagation Behavior
Lastly, Figs. 6 and 7 show how the pulse-averaged intensity of a twisted STOV beam evolves as it propagates in free space. The figures show the theoretical and simulated results, respectively, and both are organized in the same manner. Rows 1-3 show I(x, t , z) for the I, PC, and C beam cases, respectively. The I(x, t , z) in each row are encoded using the same false color scale defined by the color bars at rows' end. Columns 1-6 show the I(x, t , z) versus Fresnel number N F . Fig. 8 shows the theoretical I(x, t , z) for the GSM pulsed beam and is included for comparison. To aid the reader, row and column headings have been added to Figs. 6-8. Considering the quality of the source-plane results presented in Figs. 3-5, it is not surprising that Figs. 6 and 7 are in excellent agreement. In addition, we can see the effect of space-time coupling by comparing the I(x, t , z) in Figs. 6 and 7 to the corresponding GSM pulsed beam intensities in Fig. 8. The GSM pulsed beam, with no spatiotemporal coupling, is unchanged and diffracts in the t and x dimensions, respectively. As a result, the orientations of the Gaussian ellipses in Fig. 8 do not change. On the other hand, for the twisted STOV beam, spatiotemporal coupling causes diffraction to affect both the beam's temporal and spatial pulse widths. The relationship between these widths and the subsequent redistribution of intensity manifests as beam rotation. This can  be seen in Figs. 6 and 7, where the ellipse orientations rotate counterclockwise as N F decreases or z increases.
In all twisted STOV beam cases, the on-axis null (due to the second-order STOV) fades as a result of partial coherence combined with propagation. As physically expected, the fading is the slowest in the C beam case, where the effects of the coherent on-axis STOV are longer lasting. Indeed, starting with the C beam N F = 2.5 column and proceeding to the right, we observe two intensity minima (not complete nulls) produced by the splitting of the on-axis second-order STOV into two first-order vortices. These details are lost in the I and PC results because of lower temporal and spatial coherence.

Physical Beam Synthesis
Before concluding, we discuss how to physically generate twisted STOV beams. A twisted STOV field realization can be synthesized using a Fourier transform pulse shaper (FTPS) [60], [66]- [69]. This optical device consists of two identical diffraction gratings separated by a 4f lens system composed of two cylindrical lenses (CLs). At the center of the 4f system, i.e., in the back focal plane of CL 1 and in the front focal plane of CL 2, is a spatial light modulator (SLM).
Assuming a coherent pulsed beam is input into the FTPS, the first grating-CL combination performs a Fourier transform, mapping the input beam's spectrum into physical space at the plane of the SLM. The SLM modifies the field in the x-ω domain, producing a Fourier transformed instance of a twisted STOV field. Note that since the SLM operates in the x-ω domain, a Fourier transform must be performed on the t dimension of the kernel H in Eq. (2). Lastly, the second grating-CL combination again Fourier transforms the field, reversing the spectrum-to-space mapping of the first grating-CL system, resulting in a twisted STOV instance similar to those in Fig. 1. Partial coherence manifests by incoherently summing many such realizations, or pulses.
Ideally, the SLM would operate at the source's pulse repetition frequency (PRF). This would result in a sequence of statistically independent pulses, and therefore, quick convergence to the desired twisted STOV beam. If the SLM has a refresh rate less than the source's PRF, pulses in the sequence will be correlated. Nevertheless, the desired beam can still be realized just by integrating more pulses [69].

Conclusion
In this paper, we introduced twisted STOV beams. These beams possessed a coherent optical vortex and a random twist coupling their space and time dimensions, resulting in a stochastic pulsed beam with controllable partial coherence and transverse OAM. This latter characteristic distinguished twisted STOV fields from the more common spatial vortex and twisted beams, which possess longitudinal OAM. This difference in OAM direction creates possibilities for twisted STOV beams to be utilized in beam-control applications, such as optical tweezing, atomic optics, optical communications, and laser ablation, in novel ways.
We began the paper with the mathematical derivation of the MCF for twisted STOV beams. We then proceeded to derive relations for the linear and angular momentum densities, as well as an integral expression for the propagation of the twisted STOV MCF (ultimately evaluated numerically). We simulated the generation of twisted STOV beams and investigated their free-space propagation properties. We compared the theoretical pulse-averaged intensity and 2D cuts through the 4D MCF to the corresponding statistical moments computed from Monte Carlo simulations. The results were found to be in good agreement, proving that we had successfully generated twisted STOV fields with the desired statistics.
We lastly described how to physically realize twisted STOV beams. Constructing a Fourier transform pulse shaper, generating field realizations, and experimenting with twisted STOV beams are the next steps in this research.