Generating Electromagnetic Schell-Model Sources Using Complex Screens with Spatially Varying Auto- and Cross-Correlation Functions

We present a method to generate any physically realizable electromagnetic Schell-model source. Our technique can be directly implemented on existing vector-beam generators that utilize spatial light modulators for coherence control, beam shaping, and relative phasing. This work significantly extends published research on the subject, where control over the partially coherent source’s cross-spectral density matrix was limited. We begin by presenting the statistical optics theory necessary to derive and implement our method. We then apply our technique, both analytically and in simulation, to produce two electromagnetic Schell-model sources from the literature. We demonstrate control over the full cross-spectral density matrices of both partially coherent beams. We compare the simulated results of these two sources to the corresponding theoretical or designed quantities to validate our approach. We find, through examination of the two-dimensional correlation coefficients, that both sources converge to their desired, ensemble (or by ergodicity, “long-time”) statistics within 500 random field instances. Our method and subsequent findings will be useful in any application where control over beam shape, polarization, and spatial coherence are important. These include but are not limited to free-space/underwater optical communications, directed energy, medicine, atomic optics, and optical tweezers.

Several recent papers have presented theory and some experimental results where complex VSMSs have been produced using the above approach [3,[17][18][19][20]25,22,[26][27][28]. A key theoretical problem present in the cited works has been control over the spatial cross-correlation function of the orthogonal child beams. As will be shown, the crosscorrelation function of the child beams is predominately determined by their individual self-or auto-correlation functions. Mathematically, this affects the off-diagonal elements of the cross-spectral density (CSD) matrix and physically limits the ability to synthesize VSMSs with S 2 and S 3 Stokes vector components [29]. For instances, Refs. [3,18,[25][26][27] discuss synthesizing electromagnetic Gaussian Schell-model sources (EGSMSs). The EGSMS parameter constraints derived therein are much more restrictive than the physical source realizability conditions derived in Refs. [3,30]. More extreme examples can be found in Refs. [22,28]. In the former, the authors find only one condition under which an electromagnetc multi-Gaussian Schell-model source (EMGSMS) can be produced with off-diagonal CSD matrix elements [22]. In the latter, the inability to control the full CSD matrix limits far-zone beam control to only two Poincaré sphere parameters-the total intensity S 0 and the degree of polarization S , 0 P and the angle of polarization ψ, or S 0 and the ellipticity angle χ [28].
In this paper, we present a way to control the full CSD matrix and therefore produce any physically realizable VSMS. Our method does not require any additional optical components or hardware and can be directly implemented using the interferometer setup highlighted above and the more detailed setups described in Refs. [3,18,22].
We begin by presenting the theory underpinning our approach. Then, as examples, we generalize the work presented in Refs. [22,28], namely, to produce any physically realizable EMGSMS and an engineered VSMS that radiates a beam with a customizable far-zone polarization state, i.e., controllable S ψ , , 0 P , and χ . We demonstrate and validate our approach with Monte Carlo simulations, where we generate an EMGSMS which, until now, has been impossible to produce, and an engineered VSMS that radiates a far-zone beam with an S ψ , , 0 P , and χ that are complex grayscale images. We compare the simulated results to the EMGSMS theoretical predictions or the desired S ψ , , 0 P , and χ images, whichever is applicable. We present the correlation coefficients versus Monte Carlo trial number to show the convergence of our method. The correlation coefficient results will be useful to those who implement our technique for a specific application. Lastly, we conclude this paper with a summary of our analysis and a brief list of potential applications.

Theory
The statistical behavior of a wide-sense stationary (WSS), planar, partially coherent source is fully described by its CSD matrix [2,31], viz., where ω is the radian frequency (hereafter assumed and suppressed),= , . E α is the stochastic complex amplitude of the optical field's α polarization component.
Here, we assume that W̲ takes a Schell-model form [2,31], such that where τ α is the deterministic complex amplitude of the α field component and μ αβ is the complex cross-correlation function of the α and β field components. When it comes to physically realizing a particular W̲ , controlling μ αβ , where ≠ α β (in more direct terms, μ xy ), has been a challenge. Overcoming this challenge is the purpose of this paper. In the sections that follow, we show that μ xy can be controlled by spatially varying the cross-correlation coefficient between the white-noise, Gaussian random numbers that seed E x and E y .

Generating stochastic vector field realizations
We begin with the following model for a stochastic vector field realization:= where T α is a random complex "screen" (for the α field component) formed from circular complex Gaussian random numbers. Taking the vector auto-correlation of Eq. (3), i.e., and comparing the result to Eq. (2) yields the equality Thus, a field instance drawn from the WSS process described by W̲ can be generated by producing two correlated complex screens, T x and T y , with correlation functions given by μ αβ . The challenge of course, is generating T x and T y when, in general, ≠ ≠ μ μ μ xx yy xy . We start with = α β, as this does affect the ≠ α β case. The = α β case concerns the diagonal elements of W̲ . Since W̲ is of Schell-model form, we can generate instances of T x and T y with the desired μ xx and μ yy , respectively, by filtering delta-correlated circular complex Gaussian random numbers. For computational efficiency, the filtering operation is commonly performed in the spatial frequency domain via the convolution theorem and the fast Fourier transform (FFT). This process has been derived in the literature many times [18,22,[32][33][34][35][36]. The final result is where r α is an × N N y x grid of zero-mean, unit-variance circular complex Gaussian random numbers, i j , are discrete spatial indices, m n , are discrete spatial frequency indices, and L L , x y are the grid dimensions in meters in the x and y directions, respectively. Lastly, Φ αα is the spatial power spectrum of T α , i.e., wherê= + f x y f f x y , via the Wiener-Khinchin theorem [31,37]. Moving on to the ≠ α β case, evaluating Eq. (6) creates a T α with the desired ensemble auto-correlation μ αα . For μ xy , we need to examine the cross-correlation of T x with T y , namely, where δ x [ ] is the discrete Dirac delta function and R is the correlation coefficient between the r x and r y random numbers [3,18,22,[25][26][27][28]. Simplifying Eq. (8) with this result produces Thus, to control μ xy , xx yy xy (11) where Φ xy is the cross-power spectrum and equal to the Fourier transform of μ xy [see Eq. (7)]. With R being constant, there is little flexibility in controlling μ xy . We present an example of this later on in the paper. Next, however, we allow R to spatially vary and discuss how to generate r x and r y (ultimately, T T , x y x y x y x y x y where the superscripts "r" and "i" refer to the real and imaginary parts, respectively. The moments in Eq. (12) must be   (13) to ensure that the generated source W̲ is of Schell-model form-more accurately, that T x and T y are jointly homogeneous or equivalently that * . Substituting the above into Eq. (8) and simplifying yields This last relation leads to an expression similar to Eq. (11), but much more general: where . The next step is to generate r x and r y with cross-correlation coefficients R r,i given in Eq. (15). This can be achieved using Cholesky decomposition [22,38,39]. We note that using Cholesky decomposition to generate r x and r y for a general source given by the covariance matrix in Eq. (12) is impractical for reasonably sized grids. For example, assume that the desired T x and T y are × 512 512 grids, which is a common size for commercially available SLMs. The covariance matrix, in this case, is a staggering × 2 2 20 20 matrix, which requires a supercomputer to store and compute the Cholesky factors.
Fortunately, there is a much more efficient way to generate r x and r y , since T x and T y are jointly homogeneous. Note that in Eq. (13), the Dirac delta functions reduce the dimensions of the cross-correlation function from four to two. Thus, for any r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r where evaluation of r r r r R , , , , x x y y r i r i r , and R i at a particular m n , is implied. Calculating the Cholesky factors for Σ̲ is a relatively straightforward task, i.e., We can now generate r x and r y by first producing four × N N y x grids of zero-mean, unit-variance, independent Gaussian random numbers-r r r , , 1 Lastly, realizations of T T , x y , and E can be produced using Eqs. (6) and (3), respectively.
In the next section, we present examples showing how to apply the above analysis to generate any physically realizable EMGSMS and an engineered VSMS that radiates a far-zone beam with controllable S ψ , , 0 P , and χ (hereafter referred to as the engineered VSMS).

EMGSMS
The CSD matrix elements of an EMGSMS are [3,7] where A α and σ α are the amplitude and r.m.s width of the α field component, δ αβ is the cross-correlation width, B αβ is the complex crosscorrelation coefficient, and C 0 is a normalization factor equal to For the above source to be physically realizable, the EMGSMS parameters must satisfy [3,7] We note that when = M 1, the EMGSMS simplifies to an EGSMS, and thus, this example includes EGSMSs as a special case.
Comparing Eq. (4) to Eq. (19), we see that where C α is a complex constant, such that We can now apply the analysis in Sections "Generating stochastic vector field realizations" and "Spatially varying R" to generate EMGSMS field instances. Using Eq. (7), the spatial power spectrum Φ αα We now find the required R r,i using Eq. (15): Considering the fork inequality on the second line of Eq. (21), ⩽ ⩽ a 0 1 , where = a 0 when = B 0 xy and = a 1 when δ xy equals its maximum physically allowed value. For g, even if δ xy equals its minimum value, the numerator of g is a faster (narrower) function of f than the denominator, and therefore, ⩽ ⩽ g f 0 ( ) 1 for all f. The above conditions on a and g imply that ⩽ ⩽ R 0 1 r . Most importantly, since all δ xy satisfying the fork inequality in Eq. (21) map to an ∈ − R [ 1,1] r,i , all physically realizable EMGSMSs can be generated using the procedure and analysis presented in the previous paragraph.
To demonstrate the significance of this contribution, we briefly discuss generating EMGSMSs using the approach in Refs. [3,25,22,18,[26][27][28]. As discussed in Section "Generating stochastic vector field realizations", in past works, Eq. (11) was used to determine a spatially invariant (or constant) R to control μ xy . Reference [22] considered an EMGSMS and derived the following condition:  Thus, this prior method can only produce a small subset of physically realizable EMGSMSs. This stands in sharp contrast to the technique developed in this paper.

Engineered VSMS
Here, we design an electromagnetic Schell-model source that radiates a far-zone beam with completely controllable Poincarè sphere parameters-S ψ , , 0 P , and χ . This extends the work in Ref. [28], where using Eq. (11), only two of the parameters could be controlled simultaneously.
We begin with the random optical field given in Eq. (3). Like in Ref. [28], we assume that = = τ τ τ x y and that E is a realization of an electromagnetic quasi-homogeneous source [3,40].
Following the procedure in Ref. [28], we take the vector auto-correlation of Eq. (3), propagate the resulting W̲ to the far-zone, and then evaluate the far-zone W̲ at = = ρ ρ ρ 1 2 producing where λ is the wavelength, z is the distance to the far-zone observation point, Φ αβ is the auto-or cross-power spectrum [recall Eq. (7)], and T is the auto-correlation of τ . The Stokes parameters in terms of the CSD matrix elements are  (30) where S 0 is the total average intensity, ⩽ ⩽ 0 1 P is the degree of polarization, − < ⩽ π ψ π /2 /2 is the angle of polarization, and − ⩽ ⩽ π χ π /4 /4 is the ellipticity angle [2,3,29]. The functional dependencies of the Stokes parameters, Poincaré sphere parameters, and CSD matrix elements on ρ and z are assumed and suppressed. Substituting in Eq. (29) and simplifying produceŝ̂= We know from Eq.
With Eq. (33), we use Eq. (18) to find r r , x y , Eq. (6) to generate T T , x y , and lastly, Eq. (3) to form a vector field realization. To observe the desired Poincarè sphere parameters, we must propagate many vector field instances to the far zone and compute the Stokes parameters averaged over the ensemble of random field realizations [see Eq. (30)]. From these Stokes parameters, we then compute S ψ , , 0 P , and χ . It is instructive to explore some special polarization cases. We start with unpolarized light, where = 0 P . We see at once from Eq.   , and = ψ π 0, /2, and if so, let R r,i be any number in − [ 1,1].

Simulation
In this section, we generate an EMGSMS and an engineered VSMS with controllable far-zone S ψ , , 0 P , and χ using the expressions derived in Sections "EMGSMS" and "Engineered VSMS", respectively. Before presenting and analyzing the results, we briefly discuss the simulation setup.

Setup
For these simulations, we used = = N N 1024 x y grids with a grid spacing = Δ 15 μm. To simulate a Meadowlark Optics Model P512-635 liquid crystal SLM [41], we "windowed" the source plane field in Eq.
where = D r e c tx 7.68mm, ( ) was the rectangle function defined in Ref. [42], and τ α was The simulated wavelength was = λ 632.8nm. The complex screens, T x and T y , were generated using the procedures described in Sections "EMGSMS" and "Engineered VSMS". The simulated EMGSMS parameters are given in Table 1. The desired S ψ , , 0 P , and χ were grayscale images (photographs) shown in Fig. 4(a), (c), (e), and (g), respectively.
We generated 50,000 EMGSMS and engineered VSMS field instances and propagated each to the far zone using FFTs [32,33]. For the EMGSMS, we compared the theoretical far-zone Stokes parameters and W x x ( , 0, , 0)   [22]. For the engineered VSMS, we compared the simulated S ψ , , 0 P , and χ to their desired images.

EMGSMS
Figs. 1 and 2 show the far-zone EMGSMS Stokes parameters and W x x ( , 0, , 0) αβ 1 2 results, respectively. In Fig. 1, the left column of images shows the theoretical (superscript "thy" in the figure captions) Stokes parameters; the right column shows the simulated (superscript "sim" in the figure captions) results. The Stokes parameters are organized along the rows-S 0 in row 1 proceeding incrementally to S 3 in row 4. Row headings have been added to aid the reader. The theoretical and simulated Stokes parameters are plotted on the same color scale defined by the color bars at the end of each row. Lastly, row 5 shows the correlation coefficient C, i.e., x y x y x y (37) where N N x y is the number of pixels in an image, k is a pixel index, = i 0, 1, 2, 3, and S i is the average value of the i th Stokes vector element.
The layout of Fig. 2 mimics the × 2 2 CSD matrix W̲ . Each "element" is labeled for the reader's convenience and composed of a × 2 2 group of images. In each group, the theoretical W x x ( , 0, , 0) The results in Figs. 1 and 2 speak for themselves. Examining Fig. 1(i) shows that the simulated Stokes parameters converge to their theoretical counterparts in roughly 500 trials. To investigate this further, Fig. 3 shows S 0 and the polarization ellipses at several locations across the far-zone beam's profile. The polarization ellipse equations can be found in Ref. [3]. Fig. 3(a) and (b) show the theoretical and simulated results (after 500 trials), respectively. The agreement between Fig. 3(a) and (b) is qualitatively excellent. Fig. 4 shows the results for the engineered VSMS. The layout of the figure is similar to Fig. 1; however, here, the Poincarè sphere parameters are shown along the rows. Row headings have again been added to aid the reader, and the theoretical, or desired (left column) and simulated (right column) S ψ , , 0 P , and χ are plotted on the same color scale defined by the color bars at rows' end. Lastly, row 5 reports C [see Eq. (37)] for each of the Poincarè sphere parameters.

Engineered VSMS
Like the EMGSMS results, the agreement between theory and simulation is excellent. Fig. 4(i) shows convergence to the desired Poincarè sphere parameters occurs within approximately 500 trials. Fig. 5 shows the engineered VSMS results corresponding to those in Fig. 3. The results in Figs. 4 and 5 as well as the EMGSMS simulation results validate the electromagnetic Schell-model source synthesis method presented in this paper.

Conclusion
In this paper, we presented a method to produce any physically realizable electromagnetic Schell-model source by controlling the full CSD matrix W̲ of the partially coherent beam. In prior works focused on producing electromagnetic Schell-model sources, there was very limited control over the off-diagonal elements of W̲ , which negatively affected producing VSMSs with S 2 and S 3 Stokes vector components.
We overcame this limitation by using a spatially varying crosscorrelation coefficient R between the random numbers that seeded the field's polarization components. We showed both analytically and through examples how a spatially varying R allowed full control over W̲ . In particular, we generalized the work presented in Refs. [22,28] by demonstrating how to produce any physically realizable EMGSMS and an engineered VSMS that radiated a far-zone beam with customizable Poincaré sphere parameters.
We validated our approach with Monte Carlo simulations, where we produced the two electromagnetic Schell-model sources discussed above, both of which required full W̲ control. For the EMGSMS, we demonstrated independent control over all EMGSMS parameters. For the engineered VSMS, we demonstrated full control over the far-zone S ψ , , 0 P , and χ , which were complex grayscale images. We examined the convergence of the stochastic vector field instances to the theoretical or desired quantities by reporting the correlation coefficients for both sources versus Monte Carlo trial number. We found that both sources converged within 500 field instances. This result will be useful to those who implement our approach for a particular application.
The method we developed in this paper will be useful in SLM-based techniques to synthesize electromagnetic Schell-model sources. Our approach does not require any additional hardware and can be directly implemented on existing setups that produce vector beams [3][4][5]18,22,[43][44][45][46]. By expanding the domain of VSMSs that can be generated using existing setups, our work increases the utility of these systems for vector-beam and beam-shaping applications such as freespace optical communications, directed energy, atomic optics, particle manipulation, and optical tweezers.