Independently Controlling Stochastic Field Realization Magnitude and Phase Statistics for the Construction of Novel Partially Coherent Sources

: In this paper, we present a method to independently control the ﬁeld and irradiance statistics of a partially coherent beam. Prior techniques focus on generating optical ﬁeld realizations whose ensemble-averaged autocorrelation matches a speciﬁed second-order ﬁeld moment known as the cross-spectral density (CSD) function. Since optical ﬁeld realizations are assumed to obey Gaussian statistics, these methods do not consider the irradiance moments, as they, by the Gaussian moment theorem, are completely determined by the ﬁeld’s ﬁrst and second moments. Our work, by including control over the irradiance statistics (in addition to the CSD function), expands existing synthesis approaches and allows for the design, modeling, and simulation of new partially coherent beams, whose underlying ﬁeld realizations are not Gaussian distributed. We start with our model for a random optical ﬁeld realization and then derive expressions relating the ensemble moments of our ﬁelds to those of the desired partially coherent beam. We describe in detail how to generate random optical ﬁeld realizations with the proper statistics. We lastly generate two example partially coherent beams using our method and compare the simulated ﬁeld and irradiance moments theory to validate our technique.


Introduction
The design and physical realization of spatially partially coherent light is a highly active area of beam-control research. One of the foundational papers on the subject is "Designing genuine spatial correlations functions" by F. Gori and M. Santarsiero [1]. Cited hundreds of times since publication, the paper presents a simple recipe for creating any physically realizable partially coherent source. This recipe is expressed in the form of a superposition integral, namely, where W is the source's cross-spectral density (CSD) function [2,3], p is a positive, real function, H is the kernel, and ρ =xx +ŷy. In Equation (1), the dependence of W, p, and H on radian frequency ω is assumed and suppressed. Since p and H, for the most part, can be arbitrarily chosen, Equation (1) has been used to create all manner of random optical beams. More information on these sources can be found in the subject reviews and books in References [4][5][6][7][8][9]. The superposition rule, as Equation (1) is now known, specifies the CSD function, or second-order field moment, of a partially coherent source. It does not, however, stipulate the fourth-order field moment, more specifically, the second-order moment of irradiance or covariance of irradiance. With the exception of propagation through turbulent media, where the primary goal has been to quantify how partial coherence reduces the variance of irradiance (captured in a moment called the scintillation index) [4,[10][11][12], the covariance of irradiance has been a secondary concern to designers of random sources [6,7,[13][14][15][16][17][18][19][20][21][22][23].
The reason for this is because the underlying electric fields giving rise to the desired partially coherent source are assumed (in most cases, with very good reason) to be Gaussian distributed. According to the Gaussian moment theorem, all statistical moments greater than second order simplify to sums and products of first-and second-order moments [3,4,24]. The result being that the covariance of irradiance is completely determined by the CSD function, and the scintillation index is always equal to one. Physically, this means that the stochastic field realizations comprising the partially coherent source are fully developed speckle fields, where, in the case of uniformly correlated or Schell-model sources, the average size of a speckle is given by the width of the normalized CSD function, also known as the spatial coherence function [24].
In this paper, we present a method to independently control the field and irradiance moments of a partially coherent source. This new technique creates stochastic field instances with non-Gaussian statistics, and therefore allows for the modeling and simulation of random beam scattering and propagation through turbulent media. In addition, this approach can be used to engineer new partially coherent sources with potentially interesting and exploitable characteristics.
In the next section, we describe our random field model, which is the product of a deterministic function corresponding to the beam shape and two statistically independent random processes-one corresponding to the field's phase and the other to its irradiance. We then derive relationships between the phase and irradiance moments and the desired source's CSD function and covariance of irradiance. We end this section by presenting how to generate phase and irradiance instances, with the proper statistics, to produce a random field realization statistically consistent with the desired partially coherent source.
In Section 3, we present two examples where we apply our method to produce partially coherent sources. For the first source, we generate a Gaussian Schell-model (GSM) beam [2,4]. To validate our approach, we compare the resulting field and irradiance sample moments to theory and those of another GSM beam generated using circular complex Gaussian (CCG) random numbers [17,18,20]. In the second example, we create a new partially coherent beam, where the CSD function has the form of a Lajunen and Saastamoinen nonuniformly correlated (NUC) source [25], the irradiance is log-normally distributed, and the covariance of irradiance takes a pseudo-Schell model form [26]. We again compare the sample field and irradiance moments to theory to verify that we have indeed produced the desired source. Lastly, we conclude with a summary of our work.

Stochastic Field Model
We begin this section with our model for a random optical field realization: where τ is complex, deterministic, and physically represents the beam shape, i is a real, positive random function corresponding to the irradiance, and S is a real random function representing the phase.
Here, we assume that the random processes i and S are statistically independent. This assumption, made ultimately out of mathematical necessity, is generally not accurate. Take, for instance, the case of fully developed speckle fields. In Reference [24], Goodman shows that while the field amplitude and phase are statistically independent at any one spatial point, they are not jointly independent at two points. This means that in computing U(ρ) , we could separate the expectation involving i and S into the product of two expectations; however, the same operation would be improper when computing U(ρ 1 )U * (ρ 2 ) . In this paper, we are ultimately interested in generating partially developed speckle fields (non-Gaussian-distributed fields), and unfortunately, the joint magnitude and phase statistics of such fields are generally not known. While i and S most certainly are not independent, lacking additional information, it becomes nearly impossible to proceed without assuming so. Not accounting for the correlation between i and S surely misses physical phenomena.
The important statistical moments of Equation (2) are the mean field, the CSD function, the mean irradiance, and the covariance of irradiance. These moments, in order, are To simplify these expressions further, we let the single-position moments i(ρ) and exp[jS(ρ)] be constant. Initially, this might seem like a very restrictive stipulation. However, the function τ in the field model accounts for the effects of spatially varying i and exp(jS) . Equations (3)-(6) then become where µ i is the mean of i.

CSD Function
We now focus on the CSD function given in Equation (8). The ensemble average involving S is equivalent to the joint characteristic function (JCF) of the random variables S(ρ 1 ) and S(ρ 2 ). In general, S could be drawn from any probability distribution. If S is a delta-correlated random process, then computing the JCF is trivial. It is also unnecessary, as a delta-correlated S, physically manifests as a spatially incoherent field, and thus, exp{j[S(ρ 1 ) − S(ρ 2 )]} = δ(ρ 1 − ρ 2 ), where δ is the Dirac delta function. Clearly, this case is of little interest.
For an arbitrarily correlated random process, the only closed-form JCF known to the author is the Gaussian JCF [3,24]. If other-than-Gaussian JCFs exist, they can be used in Equation (8). The analysis would proceed in the same manner as that presented here, with the specifics dependent on the form of the JCF. We therefore assume that S is Gaussian distributed, substitute the Gaussian JCF into Equation (8), and simplify, ultimately producing where σ S is the standard deviation and R S is the normalized correlation function of the phase S, respectively. We now focus on the exponential in Equation (11), hereafter referred to as M. It is clear that M plays a role in determining the coherence of the optical field-in fact, it plays the dominant role. Its form is quite restrictive in that no matter what R S we choose (subject to the constraints on autocorrelation functions), M will always be in the form of an exponential function. Therefore, to proceed, we should choose a form for R S that allows us to easily relate phase parameters, such as phase correlation length and σ S , to those of the optical field we wish to generate. Generalizing the Gaussian forms, which are replete in the literature [3,13,14,16,[27][28][29][30][31][32], we let R S be where 0 < α ≤ 2 and χ is any real, positive function such that χ(ρ, ρ) = 0, χ(ρ 1 , ρ 2 ) > 0, and χ(ρ 1 , ρ 2 ) = χ(ρ 2 , ρ 1 ). The restrictions on χ ensure that R S is a genuine autocorrelation function.
While not yet evident, the limitation on α ensures that the function p in Equation (1) is positive. Substituting Equation (12) into Equation (11) does little to simplify M. However, we can make further progress by assuming that the field has strong phase fluctuations, such that σ S ≥ 2π [3,13,14,16,[27][28][29][30][31][32]. This large σ S means that M has significant value only when , substituting this into Equation (11), and simplifying yields It is important to note that by assuming a field with strong phase fluctuations (large σ S ), M [the exponential in Equation (13)] is a "fast function" compared to i(ρ 1 )i(ρ 2 ) and ultimately determines the coherence of the field. It is therefore safe to further approximate Equation (13) as We now have achieved the desired correspondence between the autocorrelation function of the phase R S [Equation (12)] and that of the field.
Further examination of Equation (14) reveals that the CSD function, that is, the statistics of the optical field, completely depend (with the exception of a constant) on the statistics of the phase S. This, plus the fact that Equations (9) and (10) show that the irradiance statistics depend solely on the moments of i, means that we have independent control over the field U and irradiance I statistics by controlling S and i, respectively. Although we can independently control U and I statistics, because S is Gaussian distributed and the Gaussian JCF has the form given in Equation (11), we are limited to field correlation functions of exponential type, like that in Equation (14).
Before proceeding to the irradiance moments, we briefly discuss the mean field, Equation (7). With S being Gaussian distributed, U(ρ) simplifies to where µ S is the mean of S. We stipulated above that σ S ≥ 2π. This large value of σ S effectively makes U(ρ) ≈ 0. Since µ S appears only in U(ρ) , it is a free parameter. We therefore let µ S = 0.

Covariance of Irradiance
We now turn to the irradiance moments in Equations (9) and (10), starting with B I . We first note that B I (ρ, ρ) = σ 2 I (ρ), where σ 2 I is called the scintillation index [4,12]. When dealing with speckle as opposed to irradiance fluctuations due to turbulence, σ 2 I is commonly used and is referred to as the speckle contrast [24]. Introducing σ 2 I into Equation (10) and rearranging the terms yields where b I is the normalized covariance of irradiance, that is, b I (ρ, ρ) = 1. To simplify things further, we define the zero-mean random process ζ(ρ) = i(ρ) − µ i and substitute it into Equation (16), viz., The expectation of i, that is, µ i , appears only as a multiplicative constant in I(ρ) [Equation (9)], W in Equation (14), and here in Equation (17). Therefore, like µ S , it is a free parameter and can be arbitrarily chosen subject to the constraint µ i > 0. We set µ i = 1, hereafter. This simplifies the mean irradiance in Equation (9) to Equation (17) states that any B I can be produced if we generate a sequence of zeromean random numbers ζ with σ I (ρ 1 )σ I (ρ 2 )b I (ρ 1 , ρ 2 ) as its autocorrelation function. It is easy to produce such a sequence using Gaussian random numbers. However, the difficulty here is that ζ cannot be Gaussian distributed, because that would make i Gaussian distributed. Since normal random numbers can have negative values, the Gaussian probability density function (PDF) cannot model irradiance fluctuations.
To overcome this difficulty, we use a method known by several different names in the literature-the inverse CDF (cumulative distribution function) transform method [33,34], translation [35][36][37][38], and NORTA (NORmal To Anything) [39,40]. The method starts with a correlated Gaussian random number sequence and maps it, via a non-linear transform, to a correlated non-Gaussian sequence called the target process [33][34][35][36][37][38][39][40]. If the number of discrete points in the underlying Gaussian sequence is large, translation (NORTA, etc.) yields a target realization with the desired PDF. The non-linear transform does distort the target correlation function, and there are many techniques (mostly iterative) to correct this [35][36][37][38]. Yura and Hanson [33,34] showed that target correlation function distortion caused by translation was minor. As their results attest, translation (without subsequent correction) yields highly accurate target realizations for many different PDFs and correlation functions. Here, we use Yura and Hanson's approach to generate i. The details are provided in the following section.

Generating S and i
In this section, we discuss how to generate realizations of S and i with the proper statistics. We start with phase S.
Following and modifying as necessary the analysis presented in Reference [41], an instance of S can be synthesized by evaluating the integral where "Re" denotes the real part and r S is a sample function drawn from a delta-correlated, zero-mean, CCG random process. Taking the autocorrelation of Equation (20) yields a superposition integral in the form of Equation (1), namely, Proof of this is shown in Appendix A. This is an important result as it proves that Equation (20) generates an S with the proper statistics.
Since Equation (20) must ultimately be evaluated numerically, we express it in discrete form: where is the Hadamard product, ∆v x and ∆v y are the grid spacings in the v x and v y directions, and ij and mn are discrete double indices corresponding to every combination of (x, y) and v x , v y coordinates, respectively. Equation (22) is a matrix-vector product, in which H S is a N x N y × N v x N v y matrix, where N x , N y , N v x , and N v y are the number of grid points in the dimensions denoted in the subscripts; p S and r S are N v x N v y × 1 vectors; and S is an N x N y × 1 vector. The output vector S must be reshaped to an N y × N x matrix in order to physically represent the two-dimensional phase of the field.
With S now complete, we turn our attention to generating realizations of i. As discussed at the end of the previous section, we use a process known as translation to synthesize i realizations from correlated sequences of Gaussian random numbers. These Gaussian sequences have the same first and second moments as ζ, namely, For this reason, we refer to them as ζ G for the remainder of the paper. Like the phase S, the autocorrelation function of ζ G is genuine, and therefore, a ζ G realization can be produced by evaluating where r ζ G , like r S , is a sample function drawn from a delta-correlated, zero-mean, CCG random process. The discrete representation of Equation (24) is identical to Equation (22) with appropriate substitutions. Having a ζ G realization, the next step is to form i by translating ζ G . This is accomplished using the non-linear transform where F G is the Gaussian CDF and F −1 T is the inverse target CDF [33][34][35][36][37][38][39][40]. In Equation (25), both F G and F −1 T are parameterized by their respective means and standard deviations (assuming, in the latter case, they exist). Any probability distribution i ≥ 0 can be selected as the target PDF.
The final step is to form a random field instance U by substituting the S and i realizations into Equation (2). In the next section, we use the above procedure to generate two example partially coherent beams.

Results and Discussion
Here, we use the analysis from Section 2 to generate two partially coherent sources. The first is a GSM beam, which we generate in two different ways-one using the results from Section 2 (which we call the P&I technique hereafter) and the second using an established method that filters CCG random numbers (reviewed in Appendix B and called the complex screen technique) [17,18,20]. We then compare the simulation results to GSM theory. The details are presented in the next section.
The second partially coherent source is a novel stochastic beam, where the CSD function has the form of a Lajunen and Saastamoinen NUC source [25], the irradiance is log-normally distributed, and the covariance of irradiance takes a pseudo-Schell model form [26]. For brevity, we refer to this source as the NUC beam. We compare the simulated W, B I , field PDFs, and irradiance PDF to theory to prove that we have successfully generated the desired NUC beam. The details of this simulation are presented after the GSM beam results.

GSM Beam
Before presenting the results, we provide a brief background on GSM sources and discuss details of the simulation. A GSM beam has the following CSD function W: where σ is the beam radius and δ is the correlation or coherence radius [2,4]. If the underlying optical fields are Gaussian distributed, the irradiance follows a negative exponential PDF, that is, where I = µ i = 1 [3,12,24]. The covariance of irradiance B I can easily be derived by applying the Gaussian moment theorem [3,4,24]: To form GSM field realizations using Equation (2), we compare Equation (26) to Equation (14) and conclude that Furthermore, α = 2, χ(ρ 1 , ρ 2 ) = |ρ 1 − ρ 2 |/δ S , where δ S is the phase correlation radius, and δ 2 S = 2σ 2 S δ 2 . We then use these relations to find p S and H S in Equation (20), namely, With H S being the Fourier kernel, Equation (20) equals the Fourier transform of the product of two Fourier transforms and simplifies to a convolution integral.
Similarly for i, we compare Equation (28) to Equation (23) and find that σ 2 I = 1 and The space-or shift-invariant form of b I in Equation (31) makes p ζ G and H ζ G in Equation (24) equal to As a result of H ζ G equaling the Fourier kernel, Equation (24), like Equation (20), simplifies to a convolution integral. Lastly, we use Equation (25) T is the inverse exponential CDF with σ I = 1.
For the simulations, we used grids that were N = 512 points per side. The spacings in the x-y and v x -v y domains were ∆x = ∆y = 97.66 µm and ∆v x = ∆v y = 40π m −1 , respectively. The GSM beam parameters were σ = 6.25 mm and δ = 1.67 mm. With σ S = 2π, δ S = √ 2σ S δ = 0.0148 m. We generated S, i, and complex screen (CS) T realizations by evaluating Equations (20), (24), and (A7) using fast Fourier transforms. These were then substituted into Equation (2) and Equation (A5) to form GSM P&I and CS field instances. From 10,000 GSM fields each, we computed the field PDFs, irradiance PDFs, and planar (two-dimensional) cuts through W and B I . The simulation computer code (MATLAB® version R2018b .m files) is included as Supplementary Materials to this paper. Figures 1-3 show the GSM simulation results. Figure 1a-d shows a CS and P&I GSM field realization. Clearly, CS and P&I fields are different. The way that P&I fields are formed means that both S and i are continuous, independent random functions. On the other hand, CS fields are fully developed speckle fields, and while the fields themselves are continuous, the phase, in general, is not. At locations where a CS field vanishes, a phase singularity, known as a branch point, occurs [3,42]. Since the phase S in a P&I field is defined at every location, even at locations where i = 0, branch points can only ever appear in propagated P&I fields. Although branch points do not appear in source-plane P&I fields, CS and P&I fields are quantitatively identical in terms of the key statistical moments discussed earlier in the paper. These results are presented in the remaining plots and figures: Figure 1e,f shows the optical field PDFs, and Figures 2 and 3 show the CSD function and irradiance results, respectively. In Figure 2a

NUC Beam
We now discuss the NUC beam simulation and results. The NUC beam CSD function is where γ is a real, two-dimensional vector that shifts the maximum of the correlation function away from the origin [25]. The other symbols are defined above. The PDF of irradiance is log-normal, viz., where the log-normal PDF parameters µ and σ are related to the i moments by Lastly, the covariance of irradiance takes the form where sinc(x) = sin(x)/x [26].  (2), respectively. The theoretical irradiance PDF and covariance of irradiance are given in Equations (27) and (28), respectively.
Following the same procedure used to generate a GSM field realization, we find, after comparing Equations (33) and (14), that The second line of Equation (37) implies α = 2, and δ 4 S = σ 2 S δ 4 . The p S and H S in Equation (20) are In contrast to the GSM beam, and Schell-model beams in general, H S is not the Fourier kernel, and therefore, Equation (20) must be evaluated as a matrix-vector product. This is generally a daunting task because of the size of H S . Here however, both H S and p S are one-dimensional in the v domain. This significantly reduces the size of H S and simplifies Equation (20) to For i, the autocorrelation function of ζ G is given in Equation (36). The p ζ G and H ζ G in Equation (24) are where rect(x) is the rectangle function defined in Reference [43]. Like p S and H S , p ζ G and H ζ G are one-dimensional in the v domain. Equation (24) simplifies to The last step is to transform ζ G to i using Equation (25), where F −1 T is the inverse log-normal CDF with µ and σ given in Equation (35).
For the NUC beam simulations, we used spatial domain grids that were 512 points per side with a spacing equal to 97.66 µm. To discretize v x , we used 1024 points with a spacing equal to 75.47 m −2 or m −1 for S or i, respectively. The NUC beam parameters were σ = 6.25 mm, δ = 5 mm, γ = −x4.17 −ŷ5 mm, σ I (ρ) = 0.5, and δ I = 1 mm. With σ S = 2π, δ S = √ σ S δ = 0.0125 m. Like the GSM simulations, we generated 10,000 P&I NUC beam field realizations and from these, computed the field PDFs, irradiance PDF, W, and B I . The MATLAB® scripts are included as Supplementary Materials.
The NUC beam simulation results are shown in Figures 4-6. With the exception of the CS results, the layout of the figures is identical to Figures 1-3. The agreement between the theory and P&I results is again excellent. We also note that the optical field PDFs in Figure 4c,d are clearly non-Gaussian. The exact forms of these functions are unknown; we found the theoretical PDFs numerically.

Conclusions
In this paper, we demonstrated a new technique to synthesize partially coherent sources by independently controlling the phase and magnitude statistics of the underlying optical fields. In so doing, we generated non-Gaussian random field instances, which can be used to model interesting optical phenomena, such as random beam scattering and propagation through turbulent media, and to design new partially coherent beams with potentially interesting characteristics.
We began the analysis by presenting our model for a random field realization. Our model consisted of three functions-one deterministic function corresponding to the beam shape and two statistically independent random functions representing the phase and irradiance. We then computed ensemble moments of our field model to relate the phase and irradiance random function statistics to those of the desired partially coherent beam. We ended the analysis by describing how to generate phase and irradiance realizations that, when substituted into our model, produced a random optical field with statistics corresponding to the desired partially coherent source.
To validate our method, we generated two partially coherent sources-a GSM beam and a novel source called a NUC beam. For the GSM source, we compared simulated field and irradiance moments computed from 10,000 field realizations to the theoretical expressions. In addition, we generated GSM field instances using the complex screen approach to compare our method to an established synthesis technique. While the field realizations generated using our approach differed from complex screen fields, the agreement amongst the statistical moments was excellent. For the NUC beam, we again generated 10,000 field realizations and compared the Monte Carlo moments to the desired, theoretical expressions. The agreement was again excellent.
In conclusion, the method presented in this paper expands existing partially coherent beam synthesis approaches by permitting control over the irradiance moments of random beams. In addition, this work has the potential to be used in the design of new partially coherent sources for beam-control applications.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, MATLAB® version R2018b (The MathWorks, Inc., Natick, MA, USA) scripts (.m files) used to perform the simulations described in the paper.

Conflicts of Interest:
The author declares no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript: CCG circular complex Gaussian CDF cumulative distribution function CS complex screen CSD cross-spectral density GSM Gaussian Schell-model JCF joint characteristic function NUC nonuniformly correlated NORTA NORmal To Anything PDF probability density function Appendix A. Proof of Equation (21) Starting with Equation (20), we bring the real-part operator inside the integrals to yield where the superscripts "r" and "i" denote real part and imaginary part, respectively. Taking the autocorrelation of Equation (A1) and bringing the expectations inside the integrals produces Recall that r S is a sample function drawn from a delta-correlated, zero-mean, CCG random process; thus, the above auto-and cross-correlations are This simplifies Equation (A2) to which is equivalent to Equation (21). Since Equation (24) is of the same form as Equation (20), this result holds for ζ G as well.

Appendix B. Generating GSM Beams Using the Complex Screen Technique
A GSM complex-screen (CS) field realization takes the form where T is the CS and a zero-mean, unit-variance, correlated CCG random process [20].
Taking the autocorrelation of Equation (A5) and comparing that result to the GSM CSD function given in Equation (26) yields the relation Clearly, the autocorrelation function of T is genuine, and therefore, we can use Equation (1) to generate realizations of T. Reference [41] provides the details for general sources. Specializing to GSM sources, a realization of T can be generated by evaluating where r is a sample function drawn from a delta-correlated, zero-mean, CCG random process. Equation (A7) is the Fourier transform of the product of two Fourier transforms, and therefore, simplifies to a convolution integral. Since r is Gaussian and the operation in Equation (A7) is linear, T and therefore U in Equation (A5) are Gaussian distributed. It then follows that the irradiance PDF and covariance of irradiance are given by Equations (27) and (28), respectively.