In this section, we theoretically introduce twisted space-frequency and space-time partially coherent sources.

### Twisted space-frequency partially coherent sources

Our analysis begins with the necessary and sufficient criterion for a genuine CSD function^{32,33}:

$$begin{aligned} Wleft( x_1,x_2,omega _1,omega _2right) = iint _{-infty }^{infty }pleft( v_x,v_omega right) Hleft( x_1, omega _1; v_x,v_omega right) H^*left( x_2, omega _2; v_x,v_omega right) text {d}v_xtext {d}v_omega , end{aligned}$$

(1)

where (omega ) is the radian frequency, *p* is any positive function, and *H* is an arbitrary kernel. Equation (1) is also referred to as the superposition rule in the literature. For simplicity, we restrict our analysis to one spatial dimension *x*.

Adapting Mei and Korotkova’s^{14} twist kernel, we let *H* be

$$begin{aligned} Hleft( x,omega ;v_x,v_omega right)= & {} Hleft( x,omega ;v_xright) Hleft( x,omega ;v_omega right) nonumber \ Hleft( x,omega ;v_xright)= & {} exp left( -frac{sigma _x}{2} x^2right) exp left( -frac{sigma _omega }{2} {bar{omega }}^2right) exp left[ -text {j}left( x – text {j}alpha mu omega right) v_xright] nonumber \ Hleft( x,omega ;v_omega right)= & {} exp left( -frac{sigma _x}{2} x^2right) exp left( -frac{sigma _omega }{2} {bar{omega }}^2right) exp left[ -text {j}left( omega – text {j}beta mu xright) v_omega right] , end{aligned}$$

(2)

where ({bar{omega }} = omega – omega _c) and (omega _c) is the radian frequency of the light, or carrier wave. We will discuss (sigma _x), (sigma _omega ), (alpha ), (beta ), and (mu ) later on in the paper. We note that other twist kernels exist in the literature^{10,15,18,19,21} and can be adapted in a similar manner as above to produce twisted space-frequency and space-time beams.

To generate twisted GSM beams, we choose *p* to be^{10,14,15,18,19,21}

$$begin{aligned} pleft( v_x,v_omega right) = pleft( v_xright) pleft( v_omega right) = sqrt{frac{alpha }{pi }}exp left( -alpha v_x^2right) sqrt{frac{beta }{pi }}exp left( -beta v_omega ^2right) . end{aligned}$$

(3)

Like *H* in Eq. (2), other *p* can be used, e.g., the multi-Gaussian *p* in Refs.^{14,34}. Substituting Eqs. (2) and (3) into (1) and evaluating the integrals yields a CSD of the form

$$begin{aligned} Wleft( x_1,x_2,omega _1,omega _2right)= & {} exp left( -frac{x_1^2+x_2^2}{4W_x^2}right) exp left[ -frac{left( x_1-x_2right) ^2}{2delta _x^2}right] nonumber \& times exp left( -frac{{bar{omega }}_1^2+{bar{omega }}_2^2}{4W_omega ^2}right) exp left[ -frac{left( {bar{omega }}_1-{bar{omega }}_2right) ^2}{2delta _omega ^2}right] exp left[ text {j}mu left( x_1 {bar{omega }}_2 – x_2 {bar{omega }}_1 right) right] , end{aligned}$$

(4)

where (W_x) and (W_omega ) are the beam and spectral pulse widths, (delta _x) and (delta _omega ) are the spatial coherence and spectral coherence widths, and (mu ) is the twist parameter. These beam parameters are not independent and are linked in a complex, nonlinear way. Referring back to Eq. (2),

$$begin{aligned} dfrac{1}{4W_x^2}&= sigma _x – dfrac{beta mu ^2}{2}qquad dfrac{1}{4W_omega ^2} = sigma _omega – dfrac{alpha mu ^2}{2} nonumber \ dfrac{1}{2 delta _x^2}&= dfrac{1}{4alpha } + dfrac{beta mu ^2}{4} qquad dfrac{1}{2 delta _omega ^2} = dfrac{1}{4beta } + dfrac{alpha mu ^2}{4} . end{aligned}$$

(5)

In addition, (left| mu right| delta _x delta _omega le 1)^{31,35,36}, and therefore, the space-frequency twist necessarily disappears in the coherent limit (delta _x, delta _omega rightarrow infty ).

Equation (4) has the same basic form as a spatially twisted GSM beam^{9,10,12,14,37}; however, here, space and frequency are statistically twisted. It is well known that the spectral density of a spatially twisted stochastic source rotates as the beam propagates. This rotation is in the plane orthogonal to the propagation direction, e.g., *x*–*y* plane for a *z* propagating wave. Twisted space-frequency beams also rotate—this time, in the *x*–(omega ) plane.

The paraxial, twisted space-frequency GSM CSD at any propagation plane (z > 0) can be found using the two-frequency Fresnel integral, namely,

$$begin{aligned} Wleft( x_1,x_2,omega _1,omega _2,zright)&= frac{exp left[ text {j}left( k_1-k_2right) zright] exp left[ frac{text {j}}{2z}left( k_1 x_2^2-k_2 x_2^2right) right] }{zsqrt{lambda _1 lambda _2}} nonumber \&quad times iint _{-infty }^{infty } Wleft( x_1^{‘},x_2^{‘},omega _1,omega _2right) exp left[ frac{text {j}}{2z}left( k_1 x^{‘}_{2}{^{2}}-k_2 x^{‘}_{2}{^2}right) right] nonumber \&quadtimes exp left[ -frac{text {j}}{z}left( k_1 x_1 x_1^{‘} – k_2 x_2 x_2^{‘}right) right] text {d}x_1^{‘} text {d}x_2^{‘}, end{aligned}$$

(6)

where (k_{1,2} = omega _{1,2}/c) and (lambda _{1,2} = 2pi /k_{1,2}) are the wavenumbers and wavelengths associated with (omega _1) and (omega _2), respectively and *c* is the speed of light. Substituting in Eq. (4) and after much calculus and algebra,

$$begin{aligned} Wleft( x_1,x_2,omega _1,omega _2,zright)= & {} frac{sqrt{k_1 k_2}}{2 z Omega } exp left[ text {j}left( k_1-k_2right) zright] exp left( -frac{{bar{omega }}_1^2+{bar{omega }}_2^2}{4W_omega ^2}right) exp left[ -frac{left( {bar{omega }}_1-{bar{omega }}_2right) ^2}{2delta _omega ^2}right] nonumber \× exp left[ frac{text {j} mu ^2}{8 z Omega ^2}left( k_1 {bar{omega }}_1^2 – k_2 {bar{omega }}_2^2 right) right] exp left[ frac{text {j}}{{2z}} left( 1 – frac{k_1 k_2}{4 z^2 Omega ^2}right) left( k_1 x_1^2 – k_2 x_2^2right) right] nonumber \& nonumber \× exp left[ -frac{left( mu {bar{omega }}_1 – k_2 x_2/z right) ^2 + left( mu {bar{omega }}_2 – k_1 x_1/z right) ^2}{16 W_x^2 Omega ^2} right] nonumber \& times exp left{ – frac{left[ left( mu {bar{omega }}_1 – k_2 x_2/z right) – left( mu {bar{omega }}_2 – k_1 x_1/z right) right] ^2}{8 delta _x^2 Omega ^2}right} nonumber \× exp left[ text {j} mu frac{k_1 k_2}{4 z^2 Omega ^2} left( x_1 {bar{omega }}_2 – x_2 {bar{omega }}_1 right) right] , end{aligned}$$

(7)

where (Omega ^2 = left[ sigma – text {j}k_1/left( 2zright) right] left[ sigma + text {j}k_2/left( 2zright) right] – left[ 1/left( 2delta _x^2right) right] ^2) and (sigma = 1/left( 4W_x^2right) + 1/left( 2delta _x^2right) ).

The exponentials on the third, fourth, and fifth lines of Eq. (7) correspond to the beam shape, coherence, and twist, respectively. Because of the initial space-frequency coupling, the spectral content of the beam affects its spatial distribution and equivalently, vice versa.

The spectral density *S* of the source can be found by evaluating Eq. (7) at the same space and frequency points^{30,31}, i.e.,

$$begin{aligned} Sleft( x,omega ,zright)= & {} Wleft( x,x,omega ,omega ,zright) nonumber \= & {} frac{N_F}{sqrt{1+4gamma _x^2+N_F^2}} exp left[ -left( 1 + mu ^2frac{4W_x^2 W_omega ^2}{1+4gamma _x^2+N_F^2}right) frac{{bar{omega }}^2}{2W_omega ^2}right] nonumber \& times exp left[ -left( frac{N_F^2}{1+4gamma _x^2+N_F^2}right) frac{x^2}{2W_x^2}right] exp left( mu frac{2 N_F}{1+4gamma _x^2+N_F^2}{bar{omega }}xright) , end{aligned}$$

(8)

where (gamma _x^2 = W_x^2/delta _x^2) and (N_F = 2 k W_x^2/z) is the coherent Gaussian beam Fresnel number. In order, the exponentials in Eq. (8) physically correspond to the spectral beam shape, spatial beam shape, and *x*–(omega ) plane rotation. We note that the coefficient in the spectral beam shape exponential is greater than or equal to one. This, when combined with the fact that the spatial beam shape is only affected by diffraction (depends on Fresnel number, spatial beam size, and coherence radius), means that the beam essentially “trades” spectrum to rotate. As (z rightarrow infty ) or (N_F rightarrow 0), the spectral beam radius asymptotes (the spectral beam shape does not appreciably change), diffraction dominates, and the beam no longer rotates.

Although evident from the numerous references describing rotating coherent beams^{28,37,38,39}, it is important to point out that beam rotation does not imply partial coherence. Rotation, therefore, is a characteristic of partially coherent twisted beams, not a defining characteristic.

### Twisted space-time partially coherent sources

Similar to the approach we used above to produce twisted space-frequency sources, we can also construct twisted space-time partially coherence sources. In many respects, these sources are more physically intuitive than their twisted space-frequency counterparts, as the rotation occurs in the *x*–*t* plane. Paraxially, *t* is closely related to the propagation direction (z), and therefore, these beams rotate or tumble as the beam propagates.

Like above, we begin with the superposition rule, this time for the MCF (Gamma )^{32,33}:

$$begin{aligned} Gamma left( x_1,x_2,t_1,t_2right)= & {} int _{-infty }^{infty }pleft( v_xright) Hleft( x_1,t_1;v_xright) H^*left( x_2,t_2;v_xright) text {d}v_x int _{-infty }^{infty }pleft( v_tright) Hleft( x_1,t_1;v_tright) H^*left( x_2,t_2;v_tright) text {d}v_t , end{aligned}$$

(9)

where *p* and *H* are

$$begin{aligned} pleft( v_x,v_tright)&= pleft( v_xright) pleft( v_tright) = sqrt{frac{alpha }{pi }}exp left( -alpha v_x^2right) sqrt{frac{beta }{pi }}exp left( -beta v_t^2right) nonumber \ Hleft( x,t;v_x,v_tright)&= Hleft( x,t;v_xright) Hleft( x,t;v_tright) nonumber \ Hleft( x,t;v_xright)&= exp left( -frac{sigma _x}{2} x^2right) exp left( -frac{sigma _t}{2} t^2right) exp left( -text {j}frac{omega _c}{2} tright) exp left[ -text {j}left( x – text {j}alpha mu tright) v_xright] nonumber \ Hleft( x,t;v_tright)&= exp left( -frac{sigma _x}{2} x^2right) exp left( -frac{sigma _t}{2}t^2right) exp left( -text {j}frac{omega _c}{2} tright) exp left[ -text {j}left( t – text {j}beta mu xright) v_tright] . end{aligned}$$

(10)

Substituting the above *p* and *H* into Eq. (9) and evaluating the integrals produces an MCF of the form

$$begin{aligned} Gamma left( x_1,x_2,t_1,t_2right)= & {} exp left( -frac{x_1^2+x_2^2}{4W_x^2}right) exp left[ -frac{left( x_1-x_2right) ^2}{2delta _x^2}right] exp left[ -text {j}omega _cleft( t_1-t_2right) right] nonumber \&quad times exp left( -frac{t_1^2+t_2^2}{4W_t^2} right) exp left[ -frac{left( t_1-t_2right) ^2}{2delta _t^2}right] exp left[ text {j}mu left( x_1 t_2- x_2 t_1right) right] . end{aligned}$$

(11)

The physical source parameters (W_x), (W_t), (delta _x), (delta _t), and (mu ) are related to (sigma _x), (sigma _t), (alpha ), and (beta ) in the same way as the corresponding twisted space-frequency beam parameters—see Eq. (5). As expected, when (mu = 0), Eq. (11) simplifies to a traditional GSM pulsed beam^{29,40,41,42,43,44}.

We can propagate the MCF in Eq. (11) to any plane (z > 0) using the following integral expression:

$$begin{aligned} Gamma left( x_1,x_2,t_1,t_2,zright)= & {} frac{1}{lambda _c z} iint _{-infty }^{infty } Gamma left[ x_1′,x_2′,t_1-frac{z}{c}-frac{left( x_1-x_1’right) ^2}{2cz},t_2-frac{z}{c}-frac{left( x_2-x_2’right) ^2}{2cz}right] text {d}x_1’text {d}x_2′. end{aligned}$$

(12)

This relation is accurate in the paraxial regime, and if the source is narrowband, i.e., (omega _c gg max left( 1/W_t,1/delta _tright) ). Substituting Eqs. (11) into (12) and neglecting terms greater than second order produces

$$begin{aligned} Gamma left( x_1,x_2,t_1,t_2,zright)approx & {} frac{k_c}{2zOmega } exp left[ -text {j}omega _c left( {bar{t}}_1-{bar{t}}_2right) right] exp left( -frac{{bar{t}}_1^2 + {bar{t}}_2^2}{4W_t^2}right) exp left[ -frac{left( {bar{t}}_1 – {bar{t}}_2 right) ^2}{2delta _t^2} right] nonumber \& times exp left[ frac{text {j}k_c mu ^2}{8 z Omega ^2}left( {bar{t}}_1^2-{bar{t}}_2^2right) right] exp left[ -frac{text {j}k_c^3}{8 z^3 Omega ^2}left( x_1^2-x_2^2right) right] nonumber \& times exp left[ -frac{left( mu {bar{t}}_1 – k_c x_2/z right) ^2 + left( mu {bar{t}}_2 – k_c x_1/z right) ^2}{16W_x^2Omega ^2} right] nonumber \& times exp left{ -frac{left[ left( mu {bar{t}}_1 – k_c x_2/z right) – left( mu {bar{t}}_2 – k_c x_1/z right) right] ^2}{8delta _x^2Omega ^2} right} nonumber \& times exp left[ text {j}mu frac{k_c^2}{4 z^2 Omega ^2}left( x_1 {bar{t}}_2 – x_2 {bar{t}}_1 right) right] , end{aligned}$$

(13)

where ({bar{t}}_{1,2} = t_{1,2} – z/c – x_{1,2}^2/left( 2czright) ),

$$begin{aligned} Omega ^2 = left( frac{1}{4W_x^2}right) ^2left( 1 + 4gamma _x^2 + N_F^2right) , end{aligned}$$

(14)

and (N_F = 2 k_c W_x^2/z) is the Fresnel number at the carrier frequency.

Like the spectral density above, the time-varying, ensemble-averaged intensity can be determined by evaluating the MCF at equal space and time points:

$$begin{aligned} Ileft( x,t,zright)= & {} Gamma left( x,x,t,t,zright) nonumber \= & {} frac{N_F}{sqrt{1+4gamma _x^2+N_F^2}} exp left[ -left( 1 + mu ^2frac{4W_x^2 W_t^2}{1+4gamma _x^2+N_F^2}right) frac{{bar{t}}^2}{2W_t^2}right] nonumber \& times exp left[ -left( frac{N_F^2}{1+4gamma _x^2+N_F^2}right) frac{x^2}{2W_x^2}right] exp left( mu frac{2 N_F}{1+4gamma _x^2+N_F^2}{bar{t}}xright). end{aligned}$$

(15)

The behavior of this beam in the *x*–*t* plane is the same as that described for the twisted space-frequency GSM source in the *x*–(omega ) plane.

As briefly stated above, the time *t* paraxially corresponds to the physical propagation dimension (z) (i.e., ({bar{t}} approx t – z/c)). As such, when (mu ne 0), a twisted space-time beam tumbles as it propagates, and like beams with STOVs^{23,25,26,27}, has a component of OAM in the (pm y) direction depending on the sign of (mu ). Figure 1 and the corresponding Supplementary Video V1 show this behavior for an example twisted space-time GSM partially coherent beam.