Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Tutorials Issues in Progress Current Issue All Issues Early Posting Feature Issues

Simulating random optical fields: tutorial

Milo Wilt Hyde

Open Access Open Access

Abstract

Numerous applications—including optical communications, directed energy, remote sensing, and optical tweezing—utilize the principles of statistical optics and optical coherence theory. Simulation of these phenomena is, therefore, critical in the design of new technologies for these and other such applications. For this reason, this tutorial describes how to generate random electromagnetic field instances or realizations consistent with a given or desired cross-spectral density matrix for use in wave optics simulations. This tutorial assumes that the reader has knowledge of the fundamental principles of statistical optics and optical coherence theory. An extensive reference list is provided where the necessary background information can be found. We begin this tutorial with a brief summary of the coherent-mode representation and the superposition rule of stochastic electromagnetic fields as these foundational ideas form the basis of all known synthesis techniques. We then present optical field expressions that apply these concepts before discussing proper sampling and discretization. We finally compare and contrast coherent-mode- and superposition-rule-based synthesis approaches, discussing the pros and cons of each. As an example, we simulate the synthesis and propagation of an electromagnetic partially coherent field from the literature. We compare simulated or sample statistics to theory to verify that we have successfully produced the desired field and are capturing its propagation behaviors. All computer programs, including detailed explanations of the source code, are provided with this tutorial. We conclude with a brief summary.

1. INTRODUCTION

Statistical optics and optical coherence theory, as they are commonly understood and taught, are generalizations of geometrical (ray) and wave (physical or Fourier) optics to include stochastic light sources and propagation through or scattering from random media. Applications and technologies that employ statistical optics theory are legion and include, but are not limited to, adaptive optics, astronomy, optical coherence tomography, Fourier transform spectroscopy, intensity interferometry, optical communications, remote sensing, directed energy, and optical tweezing [1–11]. As more and more technologies exploit statistical optics and optical coherence theory, accurate simulation of these phenomena becomes paramount.

One of the ways this can be achieved is to work directly with the statistical moments of the field. For simplicity, let us assume a planar, thermal, wide-sense stationary (WSS), electromagnetic source. Most, if not all, simulations involving such light sources are concerned with evaluating the following integrals:

$$\begin{split}\textbf{W}({{{\boldsymbol r}_1},{{\boldsymbol r}_2},\omega} ) &= \iiiint _{- \infty}^\infty \textbf{G}({{{\boldsymbol r}_1},\omega ;{\boldsymbol \rho}_1^\prime} ){\textbf{W}_t}({{\boldsymbol \rho}_1^\prime ,{\boldsymbol \rho}_2^\prime ,\omega} )\\&\quad\times{\textbf{G}^\dagger}({{{\boldsymbol r}_2},\omega ;{\boldsymbol \rho}_2^\prime} ){{\rm d}^2}\rho _1^\prime {{\rm d}^2}\rho _2^\prime ,\end{split}$$
where $\textbf{W}$ is the cross-spectral density matrix (CSDM) [3,9,12–16] and ${\textbf{W}_t}$ is the $2 \times 2$ transverse CSDM defined as
$${\textbf{W}_t}\big({{\boldsymbol \rho}_1^\prime ,{\boldsymbol \rho}_2^\prime ,\omega} \big) = \left\langle {\left[{\begin{array}{*{20}{c}}{{E_x}({{\boldsymbol \rho}_1^\prime ,\omega} )}\\{{E_y}({{\boldsymbol \rho}_1^\prime ,\omega} )}\end{array}} \right]{{\left[{\begin{array}{*{20}{c}}{{E_x}({{\boldsymbol \rho}_2^\prime ,\omega} )}\\{{E_y}({{\boldsymbol \rho}_2^\prime ,\omega} )}\end{array}} \right]}^\dagger}} \right\rangle .$$
$\textbf{G}$ is the dyadic Green’s function (impulse response) describing the optical system, ${\boldsymbol r} = {\boldsymbol \rho} + {\boldsymbol {\hat z}}z = {\boldsymbol {\hat x}}x + {\boldsymbol {\hat y}}y + {\boldsymbol {\hat z}}z$, $\omega$ is the radian frequency, and $\dagger$ is the conjugate transpose. Note that the term “optical system” in this context is quite broad and includes all scenarios where light’s interaction with the environment can be modeled by $\textbf{G}$—these include actual physical systems or devices (composed of lenses, mirrors, etc.) [17,18], free-space propagation [18–22], and propagation through or scattering from complex or random media [21,23–28].

Computing the fourfold superposition integral in Eq. (1) yields the exact $3 \times 3$ CSDM at the output of the optical system; however, to do so requires storing and computing a series of extremely large matrix-vector products. The matrices representing the elements of $\textbf{G}$ are 10-dimensional, with every combination of observation (${x_1},{y_1},{z_1}$ and ${x_2},{y_2},{z_2}$) and source ($x_1^\prime ,y_1^\prime $ and $x_2^\prime ,y_2^\prime $) points arranged column- and row-wise, respectively. If $\textbf{G}$ is shift-invariant, then Eq. (1) can be computed using fast Fourier transforms (FFTs) resulting in tremendous run-time and memory savings. Nevertheless, we still must store the input and output CSDMs, where each element of ${\textbf{W}_t}$ and $\textbf{W}$ represents four- and six-dimensional functions, respectively [29–31]. The use of coherent-mode or pseudo-mode expansions of ${\textbf{W}_t}$ (more on these below) can potentially make this problem more manageable [32–47]. Generally speaking, however, the dimensionality of the simulation quickly makes this approach prohibitive.

Besides using the CSDM directly, we could work with the optical field instead. In this case, we are interested in evaluating

$${\boldsymbol E}({{\boldsymbol r},\omega} ) = {\iint_{- \infty}^\infty \textbf{G}({{\boldsymbol r},\omega ;{{\boldsymbol \rho}^\prime}} ){{\boldsymbol E}_t}({{{\boldsymbol \rho}^\prime},\omega} ){{\rm d}^2}{\rho ^\prime}} ,$$
where ${{\boldsymbol E}_t} = {\boldsymbol {\hat x}}{E_x} + {\boldsymbol {\hat y}}{E_y}$. This reduces the dimensionality of the simulation by a factor of 2. However, since the input field ${{\boldsymbol E}_t}$ is random, we must generate many ${{\boldsymbol E}_t}$ realizations and subsequently evaluate Eq. (3) many times to obtain accurate estimates of $\textbf{W}$. Although this Monte Carlo approach to the problem provides an approximation to the field’s statistical moments [in contrast to Eq. (1)], the significant savings in computational resources is generally the deciding factor in its favor. In addition, having access to the fields themselves allows us to compute statistics that might otherwise be unavailable using the CSDM directly.

It is for these reasons that this paper presents a tutorial on how to generate or synthesize optical field realizations with prescribed statistical properties. In particular, this tutorial explains how to generate electromagnetic fields that are sample functions drawn from a WSS random process described by a CSDM. These fields would then be used in simulations of optical systems, propagation through random or complex media, scattering from surfaces, etc., which are described in other computational optics texts (see [48–52]). It is important to note that this is not new research. Indeed, techniques to synthesize optical fields with prescribed correlation or coherence properties can be found throughout the literature [1,4,48,53–67]. Although numerous, these references are only a small sampling of the published papers on the subject. This tutorial is written and organized in such a manner to provide a didactic review of this extensive body of research.

Lastly, it has been the author’s experience that seeing the actual field realization provides significant insight into the underlying phenomena, more so than theory alone. Therefore, this tutorial is also meant to serve as a tool for teachers of statistical optics. The techniques discussed in this paper could be used to create computer demonstrations of fundamental statistical optics concepts, such as the van Cittert–Zernike theorem, Young’s experiment, the Michelson interferometer, and the Hanbury Brown and Twiss effect. These would then augment the associated theory presented in the classic pedagogical texts (see [6,13,20]).

In the next section, we begin by summarizing the available techniques to generate electromagnetic field realizations consistent with a desired CSDM ${\textbf{W}_t}$. We present an expression for the optical field realization, discuss proper sampling or discretization, and examine the pros and cons of each approach. In Section 3, we generate and propagate (through free-space) an example electromagnetic partially coherent source using the approaches described in the prior section. We compare simulated statistics to theory to validate that we have indeed produced the desired stochastic field. Lastly, we conclude the paper with a brief summary.

For the remainder of this tutorial, we omit the “$t$ subscript” from the electric field ${\boldsymbol E}$ and $\textbf{W}$. It should be assumed (unless explicitly stated otherwise) that ${\boldsymbol E}$ and $\textbf{W}$ include only the transverse (to $z$) components of the field.

2. SIMULATION APPROACHES

To be genuine or physical, a CSDM must be square integrable, Hermitian, and nonnegative definite, such that

$$\int_{- \infty}^\infty \iiiint _{- \infty}^\infty {\left| {{W_{\alpha \beta}}({{{\boldsymbol \rho}_1},{{\boldsymbol \rho}_2},\omega} )} \right|^2}{{\rm d}^2}{\rho _1}{{\rm d}^2}{\rho _2}{\rm d}\omega \lt \infty ,$$
$$\textbf{W}({{{\boldsymbol \rho}_1},{{\boldsymbol \rho}_2},\omega} ) = {\textbf{W}^\dagger}({{{\boldsymbol \rho}_2},{{\boldsymbol \rho}_1},\omega} ),$$
$$\iiiint _{- \infty}^\infty {{\boldsymbol g}^\dagger}({{{\boldsymbol \rho}_2},\omega} )\textbf{W}({{{\boldsymbol \rho}_1},{{\boldsymbol \rho}_2},\omega} ){\boldsymbol g}({{{\boldsymbol \rho}_1},\omega} ){{\rm d}^2}{\rho _1}{{\rm d}^2}{\rho _2} \ge 0,\!$$
where $\alpha ,\beta = x,y$ and ${\boldsymbol g}$ is any well-behaved vector function [3,9,12–16,20,68,69]. Because of these properties, $\textbf{W}$ can be decomposed in two ways, known as the coherent-mode representation (or expansion) of $\textbf{W}$ and the superposition rule, respectively. These two decompositions form the basis of all current approaches to simulate random optical fields. We discuss each in turn, focusing mostly on the latter, since superposition-rule-based approaches can be used to generate electromagnetic field realizations consistent with any genuine $\textbf{W}$. Figure 1 pictorially summarizes the simulation methods discussed in this tutorial. Hereafter, we omit the $\omega$ for brevity; the dependence of all subsequent terms on $\omega$ should be assumed unless otherwise stated.
 figure: Fig. 1.

Fig. 1. Summary of approaches for simulating random optical fields. The green and red arrows show the directions of the pros (${+}$) and cons (${-}$) comparisons. For instance, the top arrow points to the coherent-mode and superposition-rule boxes. Thus, the pros and cons listed in those boxes relate coherent-mode-based and superposition-rule-based simulation approaches. Likewise, the bottom arrow points to methods 1 and 2 in both the coherent-mode and superposition-rule boxes. Therefore, the pros and cons listed under those methods relate coherent-mode method 1 to coherent-mode method 2 and superposition-rule method 1 to superposition-rule method 2.

Download Full Size | PDF

A. Coherent-Mode Representation

Using a multi-dimensional version of Mercer’s theorem, Emil Wolf derived the coherent-mode representation of a scalar WSS partially coherent field in the early 1980s [20,70,71]. It has since been generalized to non-stationary and electromagnetic random fields [3,15,72–76]. The coherent-mode expansion of a genuine $\textbf{W}$ is

$$\textbf{W}({{{\boldsymbol \rho}_1},{{\boldsymbol \rho}_2}} ) = \sum\limits_n {\lambda _n}{{\boldsymbol \psi}_n}({{{\boldsymbol \rho}_1}} ){\boldsymbol \psi}_n^\dagger ({{{\boldsymbol \rho}_2}} ),$$
where the eigenvalues ${\lambda _n} \ge 0$ for all $n$ and the vector eigenfunctions ${{\boldsymbol \psi}_n}$ are orthonormal with respect to the inner product,
$${\iint_{- \infty}^\infty {\boldsymbol \psi}_m^\dagger ({\boldsymbol \rho} ){{\boldsymbol \psi}_n}({\boldsymbol \rho} ){{\rm d}^2}\rho = {\delta _{\textit{mn}}}.}$$
Both ${\lambda _n}$ and ${{\boldsymbol \psi}_n}$ are solutions to the vector or coupled integral equations,
$${\iint_{- \infty}^\infty \textbf{W}({{{\boldsymbol \rho}_1},{{\boldsymbol \rho}_2}} ){{\boldsymbol \psi}_n}({{{\boldsymbol \rho}_2}} ){{\rm d}^2}{\rho _2} = {\lambda _n}{{\boldsymbol \psi}_n}({{{\boldsymbol \rho}_1}} ).}$$

1. Electromagnetic Field Realizations

Coherent-mode method 1: weighted coherent mode: Using the coherent modes of $\textbf{W}$, we can construct electromagnetic field realizations in two ways. The first is quite simply

$${{\boldsymbol E}_n}({\boldsymbol \rho} ) = \sqrt {{\lambda _n}} {{\boldsymbol \psi}_n}({\boldsymbol \rho} ),$$
where the field consists of a single coherent mode weighted by $\sqrt {{\lambda _n}}$. Taking the vector outer product of Eq. (10) and summing over all $n$ reproduces Eq. (7) [63,64,77–79].

Coherent-mode method 2: sum of randomly weighted coherent modes: In the second approach, the electric field is given by a coherent sum of randomly weighted modes, viz.,

$${\boldsymbol E}({\boldsymbol \rho} ) = \sum\limits_n {a_n}{{\boldsymbol \psi}_n}({\boldsymbol \rho} ),$$
where ${a_n}$ is a random complex number [20,36,70,71,73]. To determine the required statistics of ${a_n}$, we must find the moments of ${\boldsymbol E}$. Here, we are interested in simulating a thermal or pseudo-thermal light source. Therefore, ${\boldsymbol E}$ should be described by a circular complex Gaussian (CCG) probability density function, i.e., its real and imaginary parts are Gaussian distributed (or normally distributed), independent, zero mean, and have equal variances [6,8,20]. The first stipulation on ${\boldsymbol E}$ necessitates that ${a_n}$ be Gaussian distributed. For the mean of Eq. (11) to be zero, clearly $\langle {a_n^{\rm r}} \rangle = \langle {a_n^{\rm i}} \rangle = 0$, where the superscripts “${\rm r}$” and “${\rm i}$” denote the real and imaginary parts, respectively. Finally, the correlation of ${\boldsymbol E}$ (or the CSDM), viz.,
$$\langle {{\boldsymbol E}({{{\boldsymbol \rho}_1}} ){{\boldsymbol E}^\dagger}({{{\boldsymbol \rho}_2}} )} \rangle = \sum\limits_n \sum\limits_m \langle {{a_n}a_m^*} \rangle {{\boldsymbol \psi}_n}({{{\boldsymbol \rho}_1}} ){\boldsymbol \psi}_m^\dagger ({{{\boldsymbol \rho}_2}} ),$$
must equal Eq. (7). Consequently, $\langle {{a_n}a_m^*} \rangle = {\lambda _n}{\delta _{\textit{nm}}}$, or equivalently,
$$\big\langle {{{({a_n^{\rm r}} )}^2}} \big\rangle = \big\langle {{{({a_n^{\rm i}} )}^2}} \big\rangle = {\lambda _n}/2.$$
In short, generating two sequences of independent, zero-mean, ${\lambda _n}/2$-variance normal random numbers; summing the sequences with one of them multiplied by ${\rm j} = \sqrt {- 1}$; and computing the sum in Eq. (11) generates a thermal, electromagnetic field realization consistent with $\textbf{W}$.

2. Discussion

Something that has gone unstated to this point, but otherwise should be quite obvious, is that using either Eq. (10) or Eq. (11) to generate ${\boldsymbol E}$ requires that the coherent-mode representation of the relevant $\textbf{W}$ be known. Indeed, this is the primary limitation of using coherent modes to simulate random fields since finding ${\lambda _n}$ and ${{\boldsymbol \psi}_n}$ requires solving the coupled integral equations in Eq. (9). Consequently, few ${\lambda _n}$ and ${{\boldsymbol \psi}_n}$ have been found [15,72,73,80].

Bimodal expansions: The situation improves somewhat if the random source is fully polarized or scalar, such that $\textbf{W}$ simplifies to the scalar cross-spectral density (CSD) function $W$. In this case, Eq. (9) becomes a single integral equation, which, while still very difficult to solve, is easier than solving coupled integral equations. As a result, a few more scalar coherent-mode expansions are known (see [1,3,36] for summaries of these sources). In some cases, these scalar coherent modes can be used to generate thermal electromagnetic field realizations using bimodal expansions of $\textbf{W}$ [81,82]. Bimodal expansions are an option if the desired $\textbf{W}$ has an unknown electromagnetic coherent-mode representation [i.e., a solution to Eq. (9) has not been found], yet the scalar coherent-mode expansions for ${W_{\textit{xx}}}$ and ${W_{\textit{yy}}}$ are known. An important example is an electromagnetic Gaussian Schell-model (EGSM) source [3,4,9,13,83–87]. Here, we briefly summarize this simulation approach. The derivations of the expressions presented below can be found in [3,81,82]. An electromagnetic field realization formed from the bimodal expansion of  $\textbf{W}$ is

$${\boldsymbol E}({\boldsymbol \rho} ) = \hat{{\boldsymbol x}}\sum\limits_n {a_n}\psi _n^x({\boldsymbol \rho} ) + \hat{{\boldsymbol y}}\sum\limits_m {b_m}\psi _m^y({\boldsymbol \rho} ),$$
where ${a_n}$ and ${b_m}$ are Gaussian-distributed, random complex numbers exactly like those in Eqs. (11) and (12), and the functions $\psi _n^x$ and $\psi _m^y$ are the scalar coherent modes for ${W_{\textit{xx}}}$ and ${W_{\textit{yy}}}$, respectively. The means and variances of $a_n^{\rm r}$, $a_n^{\rm i}$, $b_m^{\rm r}$, and $b_m^{\rm i}$ are
$$\begin{split}\langle {a_n^{\rm r}} \rangle &= \langle {a_n^{\rm i}} \rangle = \langle {b_m^{\rm r}} \rangle = \langle {b_m^{\rm i}} \rangle = 0,\\\big\langle {{{({a_n^{\rm r}} )}^2}} \big\rangle &= \big\langle {{{({a_n^{\rm i}} )}^2}} \big\rangle = \lambda _n^x/2,\\\big\langle {{{({b_m^{\rm r}} )}^2}} \big\rangle& = \big\langle {{{({b_m^{\rm i}} )}^2}} \big\rangle = \lambda _m^y/2,\end{split}$$
where $\lambda _n^x$ and $\lambda _m^y$ are the eigenvalues for the ${W_{\textit{xx}}}$ and ${W_{\textit{yy}}}$ coherent-mode expansions. To generate ${\boldsymbol E}$ with the desired ${W_{\textit{xy}}}$ and ${W_{\textit{yx}}}$, we also need the covariances of ${a_n}$ and ${b_m}$. These are given by
$$\begin{split}\langle {a_n^{\rm r}b_m^{\rm r}} \rangle &= \langle {a_n^{\rm i}b_m^{\rm i}} \rangle = R_{\textit{nm}}^{\rm r}/2,\\\langle {a_n^{\rm i}b_m^{\rm r}} \rangle& = - \langle {a_n^{\rm r}b_m^{\rm i}} \rangle = R_{\textit{nm}}^{\rm i}/2,\end{split}$$
where $R_{\textit{nm}}^{\rm r}$ and $R_{\textit{nm}}^{\rm i}$ are the real and imaginary parts of
$${R_{\textit{nm}}} = \iiiint _{- \infty}^\infty {W_{\textit{xy}}}({{{\boldsymbol \rho}_1},{{\boldsymbol \rho}_2}} ){\left[{\psi _n^x({{{\boldsymbol \rho}_1}} )} \right]^*}\psi _m^y({{{\boldsymbol \rho}_2}} ){{\rm d}^2}{\rho _1}{{\rm d}^2}{\rho _2}.$$

Numerically finding coherent modes: Lastly, it is worth mentioning that coherent-mode expansions can be found numerically by diagonalizing the CSD function $W$ or CSDM $\textbf{W}$. To the authors’ knowledge, this has only been applied to the former [88], but there is no fundamental reason why the same approach could not be used to find the coherent modes of an electromagnetic source. In either case, we must sample and then store in memory the four-dimensional $W$ or $\textbf{W}$ to numerically compute the coherent modes. This is not desirable and, in fact, is what we are trying to avoid by generating field realizations in the first place. Even so, it may be feasible to use this approach to compute coherent-mode representations if $W$ (or $\textbf{W}$) is banded or sparse. Such $W$ describes an important type of random field known as a Schell-model source [1,3,4,6,13,20,62,89].

In summary, being discrete and orthogonal, coherent modes are an extremely powerful simulation and computational tool. Unfortunately, because so few are known, using coherent modes to simulate random optical fields—whether they be scalar or electromagnetic—is not widely applicable.

B. Superposition Rule

The primary limitation of using coherent modes in statistical optics simulations is their scarcity. In other words, given an arbitrary genuine $\textbf{W}$, it is highly unlikely that its coherent-mode expansion is known. In contrast, using the superposition rule, we can generate electromagnetic field realizations given any genuine $\textbf{W}$. We explain how in this section.

Initially derived by Gori and Santarsiero [90] for scalar sources and then shortly thereafter generalized to electromagnetic fields [91], the superposition rule states that a $\textbf{W}$ is genuine or physical if it can be expressed as

$$\textbf{W}({{{\boldsymbol \rho}_1},{{\boldsymbol \rho}_2}} ) = {\iint_{- \infty}^\infty \textbf{H}({{{\boldsymbol \rho}_1};{\boldsymbol v}} )\textbf{p}({\boldsymbol v} ){\textbf{H}^\dagger}({{{\boldsymbol \rho}_2};{\boldsymbol v},} ){{\rm d}^2}v} ,$$
where ${\boldsymbol v} = {\boldsymbol {\hat x}}{v_x} + {\boldsymbol {\hat y}}{v_y}$. The matrix $\textbf{p}$ in Eq. (18), i.e.,
$$\textbf{p}({\boldsymbol v} ) = \left[{\begin{array}{*{20}{c}}{{p_{\textit{xx}}}({\boldsymbol v} )}&\;{{p_{\textit{xy}}}({\boldsymbol v} )}\\[2pt]{p_{\textit{xy}}^*({\boldsymbol v} )}&\;{{p_{\textit{yy}}}({\boldsymbol v} )}\end{array}} \right],$$
is Hermitian and nonnegative definite, such that ${p_{\textit{xx}}}({\boldsymbol v}) \ge 0$, ${p_{\textit{yy}}}({\boldsymbol v}) \ge 0$, and $\det [{\textbf{p}({\boldsymbol v})}] \ge 0$ for all ${\boldsymbol v}$. Lastly, $\textbf{H}$ is
$$\textbf{H}({{\boldsymbol \rho};{\boldsymbol v}} ) = \left[{\begin{array}{*{20}{c}}{{H_x}({{\boldsymbol \rho};{\boldsymbol v}} )}&0\\0&{{H_y}({{\boldsymbol \rho};{\boldsymbol v}} )}\end{array}} \right],$$
where ${H_x}$ and ${H_y}$ are arbitrary kernels [92]. Since $\textbf{p}$ and $\textbf{H}$ can, for the most part, be chosen at will, Eq. (18) is incredibly powerful and has been used to develop many partially coherent sources [3,4,9,62,91,93].

1. Electromagnetic Field Realizations

Superposition-rule method 1: pseudo-modes interpretation: The superposition rule in Eq. (18) can be thought of as a continuous version of the discrete coherent-mode expansion in Eq. (7), and like coherent modes, there are two ways to construct electromagnetic field realizations using the superposition rule. In the first, analogous to Eq. (10), the field consists of a single vector “mode” with an amplitude related to the matrix $\textbf{p}$. In contrast to ${{\boldsymbol \psi}_n}$, these “modes” are not discrete, orthogonal, nor unique; thus, they are commonly referred to as pseudo-modes [40,47,94–96]. In the “pseudo-mode interpretation” of Eq. (18), ${\boldsymbol v}$ (although continuous) can be thought of as the pseudo-mode’s index, analogous to $n$ in Eqs. (7) and (10).

In a recent paper, the authors derived two sets of vector pseudo-modes, which could be used to generate electromagnetic field realizations [67]. Both were found by decomposing the matrix $\textbf{p}$: the first diagonalized (computed the eigendecomposition of) $\textbf{p}$, while the second expressed $\textbf{p}$ as the sum of diagonal and singular matrices. Here, we summarize the latter approach since it has computational advantages over the former; namely, the field expression is numerically well-behaved for all values of ${\boldsymbol v}$—more detail can be found in [67].

Let $\textbf{p}$ be the sum of diagonal and singular matrices, such that

$$\begin{split}{\textbf{p}({\boldsymbol v} )}&={ \left[{\begin{array}{*{20}{c}}{A({\boldsymbol v} )}&0\\0&{A({\boldsymbol v} )}\end{array}} \right] + \left[{\begin{array}{*{20}{c}}{B({\boldsymbol v} )}&{C({\boldsymbol v} )}\\{{C^*}({\boldsymbol v} )}&{D({\boldsymbol v} )}\end{array}} \right]}\\ &={ \left[{\begin{array}{*{20}{c}}{A({\boldsymbol v} )}&0\\0&{A({\boldsymbol v} )}\end{array}} \right] + \left[{\begin{array}{*{20}{c}}{{u_x}({\boldsymbol v} )}\\{{u_y}({\boldsymbol v} )}\end{array}} \right]{{\left[{\begin{array}{*{20}{c}}{{u_x}({\boldsymbol v} )}\\{{u_y}({\boldsymbol v} )}\end{array}} \right]}^\dagger},}\end{split}$$
where ${u_x}$ and ${u_y}$ are
$$\begin{split}{{u_x}({\boldsymbol v} )}&={ \sqrt {B({\boldsymbol v} )} \exp \left\{{{\rm j}\frac{1}{2}\arg [{C({\boldsymbol v} )} ]} \right\}},\\{{u_y}({\boldsymbol v} )}&={ \sqrt {D({\boldsymbol v} )} \exp \left\{{- {\rm j}\frac{1}{2}\arg [{C({\boldsymbol v} )} ]} \right\}.}\end{split}$$
Substituting Eq. (21) into the superposition rule in Eq. (18) and computing the matrix multiplications yields
$$\begin{split}{\textbf{W}({{{\boldsymbol \rho}_1},{{\boldsymbol \rho}_2}} )}&={ {\iint_{- \infty}^\infty A({\boldsymbol v} )\left[{\begin{array}{*{20}{c}}{{H_x}({{{\boldsymbol \rho}_1};{\boldsymbol v}} )}\\0\end{array}} \right]{{\left[{\begin{array}{*{20}{c}}{{H_x}({{{\boldsymbol \rho}_2};{\boldsymbol v}} )}\\0\end{array}} \right]}^\dagger}{{\rm d}^2}v}}\\ &+{ {\iint_{- \infty}^\infty A({\boldsymbol v} )\left[{\begin{array}{*{20}{c}}0\\{{H_y}({{{\boldsymbol \rho}_1};{\boldsymbol v}} )}\end{array}} \right]{{\left[{\begin{array}{*{20}{c}}0\\{{H_y}({{{\boldsymbol \rho}_2};{\boldsymbol v}} )}\end{array}} \right]}^\dagger}{{\rm d}^2}v}}\\&+{ {\iint_{- \infty}^\infty \left[\!{\begin{array}{*{20}{c}}{{u_x}({\boldsymbol v} ){H_x}({{{\boldsymbol \rho}_1};{\boldsymbol v}} )}\\{{u_y}({\boldsymbol v} ){H_y}({{{\boldsymbol \rho}_1};{\boldsymbol v}} )}\end{array}} \!\right]{{\left[\!{\begin{array}{*{20}{c}}{{u_x}({\boldsymbol v} ){H_x}({{{\boldsymbol \rho}_2};{\boldsymbol v}} )}\\{{u_y}({\boldsymbol v} ){H_y}({{{\boldsymbol \rho}_2};{\boldsymbol v}} )}\end{array}} \!\right]}^\dagger}{{\rm d}^2}v} .}\end{split}$$
Note that $\textbf{W}$ is now the sum of three CSDMs, each in a factorized form similar to the coherent-mode representation in Eq. (7). Therefore, we need three independent electric fields, similar in form to Eq. (10), to reproduce $\textbf{W}$; they are
$$\begin{split}{{{\boldsymbol E}_1}({{\boldsymbol \rho};{\boldsymbol v}} )}&={ {\boldsymbol {\hat x}}\sqrt {A({\boldsymbol v} )} {H_x}({{\boldsymbol \rho};{\boldsymbol v}} )},\\{{{\boldsymbol E}_2}({{\boldsymbol \rho};{\boldsymbol v}} )}&={ {\boldsymbol {\hat y}}\sqrt {A({\boldsymbol v} )} {H_y}({{\boldsymbol \rho};{\boldsymbol v}} )},\\{{{\boldsymbol E}_3}({{\boldsymbol \rho};{\boldsymbol v}} )}&={ {\boldsymbol {\hat x}}{u_x}({\boldsymbol v} ){H_x}({{\boldsymbol \rho};{\boldsymbol v}} ) + {\boldsymbol {\hat y}}{u_y}({\boldsymbol v} ){H_y}({{\boldsymbol \rho};{\boldsymbol v}} )}\\ &={ {\boldsymbol {\hat x}}\sqrt {B({\boldsymbol v} )} \exp \left\{{{\rm j}\frac{1}{2}\arg \left[{C({\boldsymbol v} )} \right]} \right\}{H_x}({{\boldsymbol \rho};{\boldsymbol v}} )}\\ &\quad+{ {\boldsymbol {\hat y}}\sqrt {D({\boldsymbol v} )} \exp \!\left\{{- {\rm j}\frac{1}{2}\arg \left[{C({\boldsymbol v} )} \right]} \right\}{H_y}({{\boldsymbol \rho};{\boldsymbol v}} ).}\end{split}$$
Computing ${{\boldsymbol E}_1}{\boldsymbol E}_1^\dagger$, ${{\boldsymbol E}_2}{\boldsymbol E}_2^\dagger$, and ${{\boldsymbol E}_3}{\boldsymbol E}_3^\dagger$ and summing over all ${\boldsymbol v}$ yields Eq. (23). The only details that remain are the functions $A,\;B,\;C$, and $D$; they are
$$\begin{split}{A({\boldsymbol v} )}&={ \frac{1}{2}\left\{{\left[{{p_{\textit{xx}}}({\boldsymbol v} ) + {p_{\textit{yy}}}({\boldsymbol v} )} \right] - \sqrt {{{\left[{{p_{\textit{xx}}}({\boldsymbol v} ) - {p_{\textit{yy}}}({\boldsymbol v} )} \right]}^2} + 4{{\left| {{p_{\textit{xy}}}({\boldsymbol v} )} \right|}^2}}} \right\}},\\{B({\boldsymbol v} )}&={ \frac{1}{2}\left\{{[{{p_{\textit{xx}}}({\boldsymbol v} ) - {p_{\textit{yy}}}({\boldsymbol v} )} ] + \sqrt {{{\left[{{p_{\textit{xx}}}({\boldsymbol v} ) - {p_{\textit{yy}}}({\boldsymbol v} )} \right]}^2} + 4{{\left| {{p_{\textit{xy}}}({\boldsymbol v} )} \right|}^2}}} \right\}},\\{C({\boldsymbol v} )}&={ {p_{\textit{xy}}}({\boldsymbol v} )},\\{D({\boldsymbol v} )}&={ \frac{1}{2}\left\{{[{{p_{\textit{yy}}}({\boldsymbol v} ) - {p_{\textit{xx}}}({\boldsymbol v} )} ] + \sqrt {{{[{{p_{\textit{xx}}}({\boldsymbol v} ) - {p_{\textit{yy}}}({\boldsymbol v} )} ]}^2} + 4{{| {{p_{\textit{xy}}}({\boldsymbol v} )} |}^2}}} \right\}.}\end{split}$$
Superposition-rule method 2: incoherent interpretation: In addition to Eq. (24), we can also generate superposition-rule field realizations in a manner similar to Eq. (11) and the bimodal expansion in Eq. (14). This approach implements the “incoherent interpretation” of the superposition rule, which was originally presented by Shirai [92]. We begin with a spatially incoherent, partially polarized source, whose polarization or coherency matrix [6,8,13,20,21,97–99] is equal to $\textbf{p}$. The ${E_x}$ and ${E_y}$ emitted by this source pass through two linear optical systems with impulse responses ${H_x}$ and ${H_y}$, respectively. The output or filtered ${E_x}$ and ${E_y}$ form the thermal electromagnetic field realization.

Let us return to the superposition rule in Eq. (18) and write the CSDM in element form, i.e.,

$${W_{\alpha \beta}}({{{\boldsymbol \rho}_1},{{\boldsymbol \rho}_2}} ) = {\iint_{- \infty}^\infty {p_{\alpha \beta}}({\boldsymbol v} ){H_\alpha}({{{\boldsymbol \rho}_1};{\boldsymbol v}} )H_\beta ^*({{{\boldsymbol \rho}_2};{\boldsymbol v}} ){{\rm d}^2}v} .$$
It is important to realize that ${W_{\textit{xx}}}$ and ${W_{\textit{yy}}}$ must both be genuine scalar CSD functions because either can be isolated from an arbitrary $\textbf{W}$ simply by passing the field through a polarizer. As a result, we can generate ${E_\alpha}$ (${E_x}$ or ${E_y}$) realizations consistent with a given ${W_{\alpha \alpha}}$ (${W_{\textit{xx}}}$ or ${W_{\textit{yy}}}$) via
$${E_\alpha}({\boldsymbol \rho} ) = {\iint_{- \infty}^\infty {r_\alpha}({\boldsymbol v} ){{\left[{\frac{1}{2}{p_{\alpha \alpha}}({\boldsymbol v} )} \right]}^{1/2}}{H_\alpha}({{\boldsymbol \rho};{\boldsymbol v}} ){{\rm d}^2}v} ,$$
where ${r_\alpha}$ is a delta-correlated, CCG-distributed, random function and physically models the spatially incoherent source described in the prior paragraph [65,66]. Multiplying ${E_\alpha}$ by its conjugate and averaging over all possible ${E_\alpha}$ reproduces ${W_{\alpha \alpha}}$ in Eq. (26). Note that Eq. (27) is the superposition-rule version of Eqs. (11) and (14) for a single component of the electric field.

Evaluating Eq. (27) generates statistically independent ${E_x}$ and ${E_y}$ with the proper self-correlations (${W_{\textit{xx}}}$ and ${W_{\textit{yy}}}$). To generate ${E_x}$ and ${E_y}$ with the proper cross correlations (${W_{\textit{xy}}}$ and ${W_{\textit{yx}}}$), we need to investigate $\langle {{E_x}({{{\boldsymbol \rho}_1}})E_y^*({{{\boldsymbol \rho}_2}})} \rangle$, which using Eq. (27) equals

$$\begin{split}&\langle {{E_x}({{{\boldsymbol \rho}_1}} )E_y^*({{{\boldsymbol \rho}_2}} )} \rangle = \iiiint _{- \infty}^\infty \langle {{r_x}({{{\boldsymbol v}_1}} )r_y^*({{{\boldsymbol v}_2}} )} \rangle\\ &\quad\times \frac{1}{2}{{\big[{{p_{\textit{xx}}}({{{\boldsymbol v}_1}} ){p_{\textit{yy}}}({{{\boldsymbol v}_2}} )} \big]}^{1/2}} { {H_x}({{{\boldsymbol \rho}_1};{{\boldsymbol v}_1}} )H_y^*\big({{{\boldsymbol \rho}_2};{{\boldsymbol v}_2}} \big){{\rm d}^2}{v_1}{{\rm d}^2}{v_2}.}\end{split}$$
Comparing this expression to Eq. (26), we see that
$$\langle {{r_x}({{{\boldsymbol v}_1}} )r_y^*({{{\boldsymbol v}_2}} )} \rangle = \frac{{2{p_{\textit{xy}}}({{{\boldsymbol v}_1}} )}}{{{{\left[{{p_{\textit{xx}}}({{{\boldsymbol v}_1}} ){p_{\textit{yy}}}({{{\boldsymbol v}_2}} )} \right]}^{1/2}}}}\delta ({{{\boldsymbol v}_1} - {{\boldsymbol v}_2}} ),$$
where $\delta$ is the Dirac delta function. Expressing ${r_x}$ and ${r_y}$ as the sum of real and imaging parts and expanding the moment in Eq. (29) produces
$$\begin{split}{{R^{\rm r}}({\boldsymbol v} )}&={ \langle {r_x^{\rm r}({\boldsymbol v} )r_y^{\rm r}({\boldsymbol v} )} \rangle = \langle {r_x^{\rm i}({\boldsymbol v} )r_y^{\rm i}({\boldsymbol v} )} \rangle = \frac{{{\rm Re}\left[{{p_{\textit{xy}}}({\boldsymbol v} )} \right]}}{{{{\left[{{p_{\textit{xx}}}({\boldsymbol v} ){p_{\textit{yy}}}({\boldsymbol v} )} \right]}^{1/2}}}}},\\{{R^{\rm i}}({\boldsymbol v} )}&={ \langle {r_x^{\rm i}({\boldsymbol v} )r_y^{\rm r}({\boldsymbol v} )} \rangle = - \langle {r_x^{\rm r}{\boldsymbol v}r_y^{\rm i}({\boldsymbol v} )} \rangle = \frac{{{\rm Im}\left[{{p_{\textit{xy}}}({\boldsymbol v} )} \right]}}{{{{\left[{{p_{\textit{xx}}}({\boldsymbol v} ){p_{\textit{yy}}}({\boldsymbol v} )} \right]}^{1/2}}}}.}\end{split}$$

In summary, we can generate thermal ${E_x}$ and ${E_y}$ realizations by evaluating

$$\begin{split}{{E_x}({\boldsymbol \rho} )}&={ {\iint_{- \infty}^\infty {r_x}({\boldsymbol v} ){{\left[{\frac{1}{2}{p_{\textit{xx}}}({\boldsymbol v} )} \right]}^{1/2}}{H_x}({{\boldsymbol \rho};{\boldsymbol v}} ){{\rm d}^2}v}},\\{{E_y}({\boldsymbol \rho} )}&={ \iint _{- \infty}^\infty {r_y}({\boldsymbol v} ){{\left[{\frac{1}{2}{p_{\textit{yy}}}({\boldsymbol v} )} \right]}^{1/2}}{H_y}({{\boldsymbol \rho};{\boldsymbol v}} ){{\rm d}^2}v,}\end{split}$$
where ${r_x}$ and ${r_y}$ are formed from zero-mean, correlated, complex Gaussian random numbers. The covariance matrix ${\boldsymbol \Sigma}$ (${\boldsymbol v}$ dependence omitted) of ${r_x}$ and ${r_y}$ is equal to
$$\begin{split}{\boldsymbol \Sigma}& = \left[{\begin{array}{*{20}{c}}{\langle {({r_x^{\rm r})}^2} \rangle}&{\langle {r_x^{\rm r}r_x^{\rm i}} \rangle}&{\langle {r_x^{\rm r}r_y^{\rm r}} \rangle}&{\langle {r_x^{\rm r}r_y^{\rm i}} \rangle}\\[3pt]{\langle {r_x^{\rm i}r_x^{\rm r}} \rangle}&{\langle {{{({r_x^{\rm i}} )}^2}} \rangle}&{\langle {r_x^{\rm i}r_y^{\rm r}} \rangle}&{\langle {r_x^{\rm i}r_y^{\rm i}} \rangle}\\[3pt]{\langle {r_y^{\rm r}r_x^{\rm r}} \rangle}&{\langle {r_y^{\rm r}r_x^{\rm i}} \rangle}&{\langle {{{({r_y^{\rm r}} )}^2}} \rangle}&{\langle {r_y^{\rm r}r_y^{\rm i}} \rangle}\\[3pt]{\langle {r_y^{\rm i}r_x^{\rm r}} \rangle}&{\langle {r_y^{\rm i}r_x^{\rm i}} \rangle}&{\langle {r_y^{\rm i}r_y^{\rm r}} \rangle}&{\langle {{{({r_y^{\rm i}} )}^2}} \rangle}\end{array}} \right] \\&= \left[{\begin{array}{*{20}{c}}1&\;\;0&\;\;{{R^{\rm r}}}&\;\;{- {R^{\rm i}}}\\0&\;\;1&\;\;{{R^{\rm i}}}&\;\;{{R^{\rm r}}}\\{{R^{\rm r}}}&\;\;{{R^{\rm i}}}&\;\;1&\;\;0\\{- {R^{\rm i}}}&\;\;{{R^{\rm r}}}&\;\;0&\;\;1\end{array}} \right].\end{split}$$
Finally, ${r_x}$ and ${r_y}$ can be directly generated from the Cholesky factor [100,101] of ${\boldsymbol \Sigma}$, such that
$$\begin{split}{{r_x}({\boldsymbol v} )}&={ {r_1}({\boldsymbol v} ) + {\rm j}{r_2}({\boldsymbol v} )},\\{{r_y}({\boldsymbol v} )}&= {r_x}({\boldsymbol v} )\left[{{R^{\rm r}}({\boldsymbol v} ) - {\rm j}{R^{\rm i}}({\boldsymbol v} )} \right] + [{{r_3}({\boldsymbol v} ) + {\rm j}{r_4}({\boldsymbol v} )} ]\\&\quad\times\sqrt {1 - {{[{{R^{\rm r}}({\boldsymbol v} )} ]}^2} - {{[{{R^{\rm i}}({\boldsymbol v} )} ]}^2}} ,\end{split}$$
where ${r_1}$, ${r_2}$, ${r_3}$, and ${r_4}$ are statistically independent, delta-correlated, standard-normal random numbers.

2. Sampling

The superposition-rule field realizations in Eqs. (24) and (31) are continuous functions of ${\boldsymbol v}$ and ${\boldsymbol \rho}$. To use either in a simulation, we must properly sample the ${\boldsymbol v}$ and ${\boldsymbol \rho}$ domains, which entails finding the maximum spacings (${\Delta _v}$ and $\Delta$) and minimum points ($N$ and $P$) of the discrete grids representing those domains. The Shannon–Nyquist sampling criterion [48,49,102], viz.,

$$\Delta \le \frac{1}{{2{f_{{\rm max}}}}},$$
where ${f_{{\rm max}}}$ is the largest frequency component of ${\boldsymbol E}$, forms the basis of this analysis. Violating Eq. (34) will alias ${\boldsymbol E}$ and produce non-physical simulation results.

Spatial domain sampling $\Delta$ and $\boldsymbol{P}$: Let us begin with the spatial ${\boldsymbol \rho}$ domain sampling analysis. Note that only ${H_\alpha}$ in Eqs. (24) and (31) depends on ${\boldsymbol \rho}$. Determining the maximum frequency component of ${\boldsymbol E}$ in the ${\boldsymbol \rho}$ domain (${f_{{\rm max}}}$) requires that we first Fourier transform ${H_\alpha}$ (with respect to ${\boldsymbol \rho}$), i.e.,

$${\tilde H_\alpha}({{\boldsymbol f};{\boldsymbol v}} ) = \iint _{- \infty}^\infty {H_\alpha}({{\boldsymbol \rho};{\boldsymbol v}} )\exp ({- {\rm j}2\pi {\boldsymbol f} \cdot {\boldsymbol \rho}} ){{\rm d}^2}\rho .$$
The maximum radius of $| {{{\tilde H}_\alpha}({{\boldsymbol f};{\boldsymbol v}})} |$ is ${f_{{\rm max}}}$. Since the radius of $| {{{\tilde H}_\alpha}({{\boldsymbol f};{\boldsymbol v}})} |$ generally increases with $v$, we assume that ${f_{{\rm max}}}$ is located at the edge of the ${\boldsymbol v}$ grid, viz., ${\boldsymbol v} = ({{\boldsymbol {\hat x}} + {\boldsymbol {\hat y}}}){D_p}/2$, where ${D_p}$ is the maximum diameter of $| {{p_{\alpha \beta}}} |$. We define ${f_{{\rm max}}}$ and ${D_p}/2$ as the ${\boldsymbol f}$ and ${\boldsymbol v}$ locations where $| {{{\tilde H}_\alpha}} |$ and $| {{p_{\alpha \beta}}} |$ fall (and remain) below a user-defined value $\epsilon$ (e.g., 0.01 or 0.001). Using a smaller value for $\epsilon$ increases both $P$ and $N$ (discussed below), resulting in a more finely sampled ${\boldsymbol E}$. Of course, this improvement in resolution comes at the expense of computational resources and run-time. With ${f_{{\rm max}}}$, the ${\boldsymbol \rho}$-domain grid spacing $\Delta$ can be found by applying the Shannon–Nyquist criterion [see Eq. (34)]. Finally, the minimum number of ${\boldsymbol \rho}$-grid points per side is equal to
$$P \ge {\rm ceil}\left({L/\Delta} \right) = {\rm ceil}\left({2L{f_{{\rm max}}}} \right),$$
where ${\rm ceil}$ rounds $x$ up to the nearest integer and $L$ is the physical side-length of the grid, which must be large enough to fit or support ${\boldsymbol E}$. We can estimate the field’s size from the square root of the spectral density, i.e., $\sqrt {{\rm tr}[{\textbf{W}({{\boldsymbol \rho},{\boldsymbol \rho}})}]}$, where ${\rm tr}$ is the trace of the matrix.

${\boldsymbol v}$-domain sampling ${\Delta _v}$ and $\boldsymbol{N}$: Proceeding to the ${\boldsymbol v}$ domain, we follow the same procedure to find the grid spacing ${\Delta _v}$ and number of points $N$. The ${\boldsymbol E}$ in Eqs. (24) and (31) are generally the products of ${p_{\alpha \beta}}$ and ${H_\alpha}$; both are functions of ${\boldsymbol v}$. Therefore, to determine ${f_{v,{\rm max}}}$, we must Fourier transform (with respect to ${\boldsymbol v}$) the product of ${p_{\alpha \beta}}$ and ${H_\alpha}$ and find the radius of the resulting function [102–105]. Since the Fourier transform of the product of two functions is the convolution of their Fourier transforms, ${f_{v,{\rm max}}}$ is equal to the sum of the individual max frequencies of ${p_{\alpha \beta}}$ and ${H_\alpha}$, namely,

$${f_{v,{\rm max}}} = f_{v,{\rm max}}^{\,\textbf{p}} + f_{v,{\rm max}}^{\,\textbf{H}},$$
where $f_{v,{\rm max}}^{\,\textbf{p}}$ and $f_{v,{\rm max}}^\textbf{H}$ are the max radii of $| {{{\tilde p}_{\alpha \beta}}({{{\boldsymbol f}_v}})} |$ and $| {{{\tilde H}_\alpha}({{\boldsymbol \rho};{{\boldsymbol f}_v}})} |$, respectively. Since the radius of $| {{{\tilde H}_\alpha}({{\boldsymbol \rho};{{\boldsymbol f}_v}})} |$ generally increases with $\rho$, we assume, like above, that $f_{v,{\rm max}}^\textbf{H}$ is located at the edge of the ${\boldsymbol \rho}$ grid, i.e., ${\boldsymbol \rho} = ({{\boldsymbol {\hat x}} + {\boldsymbol {\hat y}}})L/2$. We define $f_{v,{\rm max}}^\textbf{p}$ and $f_{v,{\rm max}}^\textbf{H}$ using the same criterion described in the prior paragraph. After finding ${f_{v,{\rm max}}}$, the grid spacing ${\Delta _v}$ [via Eq. (34)] and minimum number of grid points per side,
$$N \ge {\rm ceil}({{D_p}/{\Delta _v}} ) = {\rm ceil}({2{D_p}{f_{v,{\rm max}}}} ),$$
follow immediately.

FFT field synthesis and sampling: It is important to mention that, for some random electromagnetic fields, ${H_\alpha}$ contains a Fourier-like kernel function. A very important example of such a field is an electromagnetic Schell-model source [1,3,4,9,13,62,90,91]. In these instances, the integrals in Eq. (31) become Fourier transforms and can be computed using FFTs. The FFT links the ${\boldsymbol v}$ and ${\boldsymbol \rho}$ grids, such that $N = P$ and

$$\begin{split}{{\Delta _v}}&={ \frac{{2\pi}}{L} = \frac{{2\pi}}{{P\Delta}} = \frac{{2\pi}}{{N\Delta}}},\\\Delta &={ \frac{{2\pi}}{{{L_v}}} = \frac{{2\pi}}{{N{\Delta _v}}} = \frac{{2\pi}}{{P{\Delta _v}}},}\end{split}$$
where ${L_v} \ge {D_p}$ is the side-length of the ${\boldsymbol v}$ grid. A relatively simple way to accommodate the FFT’s constraints into the sampling procedure described above is to first find ${\Delta _v}$ and $\Delta$ by enforcing the Shannon–Nyquist criterion in the ${\boldsymbol v}$ and ${\boldsymbol \rho}$ domains, respectively. Using Eq. (39), we then compute the ${\boldsymbol v}$- and ${\boldsymbol \rho}$-grid side-lengths, ${L_v}$ and $L$. Lastly, we find the number of grid points by $N = P = {\rm ceil}({{L_v}/{\Delta _v}}) = {\rm ceil}({L/\Delta})$. Because $N$ is rounded up by the ceiling function, we need to adjust the grid spacings (${\Delta _v} = {L_v}/N$ and $\Delta = L/N$) to satisfy Eq. (39).

Impulse response $\textbf{G}$ sampling: Following the steps listed in the prior three paragraphs will produce accurate discrete representations of the ${\boldsymbol E}$ in Eqs. (24) and (31) for any ${p_{\alpha \beta}}$ and ${H_\alpha}$. However, the purpose of most simulations is not just to synthesize the partially coherent field, but the purpose is to find ${\boldsymbol E}$ at the output of some optical system—put simply, evaluate Eq. (3). To do that accurately, we must include the Green’s function $\textbf{G}$ in the sampling analysis. It is clear from the ${\boldsymbol E}$ in Eqs. (24) and (31) that $\textbf{G}$ does not affect $N$; it does impact the number of spatial grid points $P$, however. Although the form of the output field in Eq. (3) is different depending on which ${\boldsymbol E}$ we use [Eqs. (24) or (31)], both require us to evaluate

$$\iint _{- \infty}^\infty \textbf{G}({{\boldsymbol r};{{\boldsymbol \rho}^\prime}} ){\boldsymbol H}({{{\boldsymbol \rho}^\prime};{\boldsymbol v}} ){{\rm d}^2}{\rho ^\prime},$$
where ${\boldsymbol H} = {\boldsymbol {\hat x}}{H_x} + {\boldsymbol {\hat y}}{H_y}$. Therefore, similar to the ${\boldsymbol v}$-domain sampling analysis, we need to find ${f_{{\rm max}}}$ and apply the Shannon–Nyquist criterion to the product of $\textbf{G}$ and ${\boldsymbol H}$ with respect to the source coordinates ${x^\prime},{y^\prime}$ [102–105]. The convolution theorem allows us to consider $\textbf{G}$ and ${\boldsymbol H}$ separately, such that
$${f_{{\rm max}}} = f_{{\rm max}}^{\,\textbf{G}} + f_{{\rm max}}^{\boldsymbol H},$$
where $f_{{\rm max}}^\textbf{G}$ and $f_{{\rm max}}^{\boldsymbol H}$ are the max radii of $| {{\boldsymbol {\tilde G}}({{\boldsymbol r};{\boldsymbol f}})} |$ and $| {{{\tilde H}_\alpha}({{\boldsymbol f};{\boldsymbol v}})} |$, respectively. The former is located at the edge of the output grid, i.e., ${\boldsymbol r} = ({{\boldsymbol {\hat x}} + {\boldsymbol {\hat y}} + {\boldsymbol {\hat z}}}){L_{{\rm out}}}/2$, whereas the latter is at edge of the ${\boldsymbol v}$ grid [${\boldsymbol v} = ({{\boldsymbol {\hat x}} + {\boldsymbol {\hat y}}}){D_p}/2$] as before. After finding ${f_{{\rm max}}}$, both $\Delta$ and $P$ can be computed from Eqs. (34) and (36), respectively.

The last thing we must consider is the sampling of the output grid. Ultimately, this depends on how we evaluate Eq. (40). If we elect to compute Eq. (40) directly as a matrix-vector product, the input and output grids are independent. In this case, we find ${f_{{\rm max}}}$ and apply the Shannon–Nyquist criterion with respect to the observation coordinates $x,y,z$. The procedure is identical to that described above. On the other hand, if $\textbf{G}$ is shift-invariant (e.g., free-space propagation), then it is computationally advantageous to evaluate Eq. (40) using FFTs. The FFT requires that the input and output grids have the same number of points and relates $\Delta$ and $L$ to the output-grid size ${L_{{\rm out}}}$ and spacing ${\Delta _{{\rm out}}}$, respectively [see Eq. (39)]. Consequently, there is a trade-off between the sampling in the input and output grids. Many researchers have discussed this topic, especially for simulating free-space propagation using FFTs [18,48,49,106–108]. We apply this analysis in the simulation in Section 3.

3. Discussion

As mentioned above, when using the superposition rule to generate electromagnetic field realizations, we have two choices for ${\boldsymbol E}$ given in Eqs. (24) and (31), respectively; both have pros and cons. The primary benefit of using Eq. (24) in a simulation is that Eq. (3) needs to be evaluated, at most, ${N^2}$ times. Recall that $N$ is the number of ${\boldsymbol v}$-domain grid points per side; therefore, ${N^2}$ generally equals the number of pseudo-modes. However, Eq. (24) is not a realization of a thermal electromagnetic field, and although using it in simulations will yield accurate second-order field moments, all other statistics will be non-physical. Equation (31), on the other hand, is a thermal electromagnetic field realization; all statistics computed using Eq. (31) will be physical. However, this realism comes at a price because computing those statistics requires that we repeatedly evaluate Eq. (3) using an independent realization of Eq. (31) each time, also known as a trial.

Number of trials T and residual error: We can estimate the minimum number of trials $T$ by finding the trial number at which the variance of the sample second-order statistic falls below a certain value. Since ${\boldsymbol E}$ is Gaussian distributed, $T$ can be found analytically and is actually captured in a physical statistic known as the speckle contrast [8] or scintillation index [3,27,28], depending on the context. The speckle contrast is a measure of the strength of irradiance fluctuations relative to the average irradiance, i.e., ${\cal C} = {\sigma _I}/\langle I \rangle$, where ${\sigma _I}$ is the standard deviation of the irradiance [8]. Derived by Goodman [8], the speckle contrast ${\cal C}$ after integrating $T$ independent speckle patterns is

$${\cal C} = \frac{1}{{\sqrt T}}.$$
Let us assume that, in our simulation, we would like to reduce the “noise” or randomness in the second-order field moment by 99%. Inverting Eq. (42) with ${\cal C} = 0.01$ reveals that $T \ge 10{,}000$.

Coherent modes versus pseudo-modes: In contrast to coherent modes, pseudo-modes are not unique, discrete, or orthogonal. As a result, typically many more pseudo-modes are required than coherent modes to accurately represent a random electromagnetic field. This point is best illustrated through an example: Let us assume that we are simulating a scalar Gaussian Schell-model (GSM) source. A GSM source has a known coherent-mode representation (see [3,20,32,36,109,110]) and a known $p$ and $H$ (see [4,62,64,90]). In addition, expressions for the required number of coherent modes (known as Starikov’s criterion) and pseudo-modes can be found in closed form [36,64,111]. To compute these numbers, we need to specify the r.m.s. width of the GSM source $\sigma$, its spatial coherence radius $\delta$, the side-length of our grid $L$, and finally, $\epsilon$. For $\sigma = 1\;{\rm cm} $, $\delta = \sigma /5$, $L = 10\sigma$, and $\epsilon = 0.001$, Starikov’s criterion estimates that ${11^2}$ coherent modes are required to generate the GSM source; on the other hand, ${68^2}$ pseudo-modes are required to do the same.

Evidently, coherent modes are many times more efficient (computationally) than pseudo-modes; however, there are few $W$—and even fewer $\textbf{W}$—with known coherent-mode expansions (see [1,3,36] for more information on these sources). This stands in striking contrast to pseudo-modes or more generally the superposition rule, where the $\textbf{p}$ and $\textbf{H}$ for all $\textbf{W}$ are known. For instance, letting ${H_\alpha}({{\boldsymbol \rho};{\boldsymbol v}}) = {\tau _\alpha}({\boldsymbol \rho})\exp ({{\rm j}{\boldsymbol v} \cdot {\boldsymbol \rho}})$, where ${\tau _\alpha}$ is the deterministic complex shape of the $\alpha$ field component, will produce any electromagnetic Schell-model source—its $\textbf{W}$ depending on $\textbf{p}$ [91]. The fact that all electromagnetic Schell-model sources can be generated with the same ${H_\alpha}$ clearly demonstrates the power of pseudo-modes and the superposition rule. This is simply not the case for coherent modes.

3. SIMULATION EXAMPLE

In this section, we generate thermal field realizations consistent with an electromagnetic Gaussian pseudo-Schell model (EGPSM) source [66,112–117]. Pseudo-Schell model beams belong to a broad class of partially coherent fields known as nonuniformly correlated (NUC) sources. They, along with several other NUC beams, “self-focus,” i.e., their peak irradiance or spectral density occurs at a near-field distance from the source. In addition, EGPSM beams possess an on-axis null and can be constructed with spatially varying ($x,y$ varying) polarization states that are invariant with respect to propagation distance $z$ [116,117]. These characteristics make EGPSM beams potentially useful in beam control applications such as optical tweezing and particle manipulation and, therefore, a germane random field to simulate.

A. Theoretical CSDM for an EGPSM Beam

The CSDM elements for an EGPSM beam take the following form [116,117]:

$$\begin{split}{W_{\alpha \beta}}({{{\boldsymbol \rho}_1},{{\boldsymbol \rho}_2}} ) &= {A_\alpha}\frac{{{\rho _1}}}{{{\sigma _\alpha}}}\cos ({{\phi _1} - {\theta _\alpha}} )\exp \left({- \frac{{\rho _1^2}}{{\sigma _\alpha ^2}}} \right)\\& \times{A_\beta}\frac{{{\rho _2}}}{{{\sigma _\beta}}}\cos ({{\phi _2} - {\theta _\beta}} )\exp \left({- \frac{{\rho _2^2}}{{\sigma _\beta ^2}}} \right){B_{\alpha \beta}}\\&\times\exp \left[{- \frac{{{{\left({{\rho _1} - {\rho _2}} \right)}^2}}}{{\delta _{\alpha \beta}^2}}} \right],\end{split}$$
where ${A_\alpha}$, ${\sigma _\alpha}$, and ${\theta _\alpha}$ are the amplitude, radius, and shift angle of the $\alpha$ field component; ${B_{\alpha \beta}}$ is the complex correlation coefficient between the $\alpha$ and $\beta$ field components; and lastly, ${\delta _{\alpha \beta}}$ is the radius of the cross correlation function. For the CSDM in Eq. (43) to be physical or genuine, ${B_{\alpha \beta}}$ and ${\delta _{\alpha \beta}}$ must satisfy
$$\begin{split}{B_{\textit{xx}}} &= {B_{\textit{yy}}} = 1, {B_{\textit{xy}}} = B_{\textit{yx}}^*, \left| {{B_{\textit{xy}}}} \right| \le 1, {\delta _{\textit{xy}}} = {\delta _{\textit{yx}}}, \\&\qquad\sqrt {\frac{{\delta _{\textit{xx}}^2 + \delta _{\textit{yy}}^2}}{2}} \,\le\, {\delta _{\textit{xy}}} \,\le\, \sqrt {\frac{{{\delta _{\textit{xx}}}{\delta _{\textit{yy}}}}}{{| {{B_{\textit{xy}}}} |}}} .\end{split}$$
The coherent-mode representation for an EGPSM source is unknown; however, the superposition rule $\textbf{p}$ and $\textbf{H}$ are
$$\begin{split}{{p_{\alpha \beta}}({\boldsymbol v} )} &= {{B_{\alpha \beta}}\frac{{\delta _{\alpha \beta}^2}}{{4\pi}}\exp \left({- \frac{{\delta _{\alpha \beta}^2}}{4}{v^2}} \right)},\\{{H_\alpha}({{\boldsymbol \rho};{\boldsymbol v}} )}&= {{\tau _\alpha}({\boldsymbol \rho} )h({{\boldsymbol \rho};{\boldsymbol v}} )}\\&={ {A_\alpha}\frac{\rho}{{{\sigma _\alpha}}}\cos ({\phi - {\theta _\alpha}} )\exp \!\left({- \frac{{{\rho ^2}}}{{\sigma _\alpha ^2}}} \right)\exp \left({{\rm j}\rho {v_x}} \right).}\end{split}$$
Consequently, we can use either Eq. (24) or Eq. (31) to generate field realizations. Since in this simulation we are interested in synthesizing thermal ${\boldsymbol E}$, we choose the latter.

Before we can proceed to the sampling analysis, we need to specify values for the EGPSM parameters in Eq. (43). Here, we generate field realizations of an EGPSM source with ${A_x} = 1$, ${\sigma _x} = 1\;{\rm cm} $, ${\delta _{\textit{xx}}} = {\sigma _x}/3$, ${\theta _x} = \pi /3$, ${A_y} = 1.25$, ${\sigma _y} = 1.25\;{\rm cm} $, ${\delta _{\textit{yy}}} = {\sigma _y}/5$, ${\theta _y} = - \pi /4$, ${B_{\textit{xy}}} = 0.35\exp ({- {\rm j}\pi /6})$, and ${\delta _{\textit{xy}}} = 0.45\;{\rm cm} $.

B. Sampling

1. Spatial Domain Sampling $\Delta$ and P

With ${p_{\alpha \beta}}$, ${H_\alpha}$, and specific EGPSM parameter values, we can now apply Section 2.B.2 to determine the grid spacings and minimum numbers of points per side in the ${\boldsymbol v}$- and ${\boldsymbol \rho}$-domains—${\Delta _v}$, $\Delta$, $N$, and $P$, respectively. As we did Section 2.B.2, let us start with the ${\boldsymbol \rho}$ domain. Since we are propagating the field realizations in this simulation, we need to account for $\textbf{G}$ when finding $\Delta$ and $P$. To propagate ${\boldsymbol E}$, we evaluate the Fourier transform form of the Fresnel diffraction integral [18], namely,

$$\begin{split}{E_\alpha}({{\boldsymbol \rho},z} ) &= \frac{{\exp ({{\rm j}kz} )}}{{{\rm j}\lambda z}}\exp \left({\frac{{{\rm j}k}}{{2z}}{\rho ^2}} \right) \iint_{- \infty}^\infty {E_\alpha}({{{\boldsymbol \rho}^\prime},0} )\\[-3pt]&\times\exp \left({\frac{{{\rm j}k}}{{2z}}{\rho ^{^\prime 2}}} \right)\exp \left({- {\rm j}\frac{k}{z}{\boldsymbol \rho} \cdot {{\boldsymbol \rho}^\prime}} \right){{\rm d}^2}{\rho ^\prime},\end{split}$$
where $\lambda = 1\,\,\unicode{x00B5}{\rm m}$ is the wavelength and $k = 2\pi /\lambda$, using FFTs [48,49,106–108]. The relevant parts of $\textbf{G} = ({{\boldsymbol {\hat x\hat x}} + {\boldsymbol {\hat y\hat y}}})G$ in Eq. (46) are the inner and outer quadratic phase terms or spatial chirps. These functions are not bandlimited, so the best we can do is ensure that they are properly sampled over user-defined input- (source) and output- (observation) plane regions of interest—${D_{{\rm in}}}$ and ${D_{{\rm out}}}$, respectively. This analysis has been performed elsewhere [48,49,106–108]; the result is
$$\begin{split}\Delta &\le{ \frac{{\lambda z}}{{{D_{{\rm in}}} + 2\lambda z{f_{{\rm max}}}}}},\\[-3pt]N&\ge{ \frac{{{D_{{\rm in}}} + {D_{{\rm out}}}}}{\Delta}.}\end{split}$$
Here, we set ${D_{{\rm in}}} = 10\max ({{\sigma _\alpha}})$ and ${D_{{\rm out}}} = 2{D_{{\rm in}}}$, which are easily large enough to support ${\boldsymbol E}$ in both the input and output planes.

The last thing we need to determine is ${f_{{\rm max}}}$: We Fourier transform ${H_\alpha}$ with respect to ${\boldsymbol \rho}$ and find the maximum radius of $| {{{\tilde H}_\alpha}({{\boldsymbol f};{\boldsymbol v}})} |$ evaluated at ${\boldsymbol v} = ({{\boldsymbol {\hat x}} + {\boldsymbol {\hat y}}}){D_p}/2$, where ${D_p} = 4\sqrt {- \ln \epsilon} /\min ({{\delta _{\alpha \beta}}})$. Unfortunately, ${\tilde H_\alpha}$ cannot be evaluated in closed form, and we must find ${f_{{\rm max}}}$ numerically. Using the EGPSM parameters listed above, the $f$ location at which the normalized $| {{{\tilde H}_\alpha}} |$ falls and remains below $\epsilon = 0.001$ is ${f_{{\rm max}}} \approx 680.14\;{{\rm m}^{- 1}}$. Substituting everything into Eq. (47) produces $P = 1944$ and $\Delta = 64.3\,\,\unicode{x00B5}{\rm m}$ with $z = 32.7\;{\rm m} $, which corresponds to a Fresnel number ${N_F} = {\max}^2 ({{\sigma _\alpha}})k/({2z}) \approx 15$. Note that this $P$ and $\Delta$ are sufficient for all ${N_F} \lt 15$.

 figure: Fig. 2.

Fig. 2. EGPSM source-plane Stokes parameter and speckle contrast ${\cal C}$ results: (a) and (b) ${S_0}$ theory and simulation, (c) and (d) ${S_1}$ theory and simulation, (e) and (f) ${S_2}$ theory and simulation, (g) and (h) ${S_3}$ theory and simulation, and (i) and (j) ${\cal C}$ theory and simulation.

Download Full Size | PDF

2. ${\boldsymbol v}$-Domain Sampling ${\Delta _{\textit{v}}}$and N

Let us now proceed to the $\boldsymbol{v}$-domain sampling. Recall that $\textbf{G}$ has no effect on ${\Delta _v}$ and $N$; therefore, finding ${f_{v,{\rm max}}}$ via Eq. (37) is our only concern. Both $f_{v,{\rm max}}^\textbf{p}$ and $f_{v,{\rm max}}^\textbf{H}$ can be found in closed form. The resulting expressions for ${\Delta _v}$ and $N$ are

$$\begin{split}{{\Delta _v}}&\le{ \frac{{2\pi}}{{2\max ({{\delta _{\alpha \beta}}} )\sqrt {- \ln \epsilon} + {D_{{\rm in}}}}}},\\N&\ge{ {\rm ceil}\left\{{\frac{4}{\pi}({- \ln \epsilon} )\frac{{\max ({{\delta _{\alpha \beta}}} )}}{{\min ({{\delta _{\alpha \beta}}} )}}\left[\!{1 + \frac{{{D_{{\rm in}}}}}{{2\max ({{\delta _{\alpha \beta}}} )\sqrt {- \ln \epsilon}}}} \!\right]} \right\}\!.}\end{split}$$
Substituting in the relevant EGPSM parameter values, ${D_{{\rm in}}} = 10\max ({{\sigma _\alpha}})$, and $\epsilon = 0.001$ yields $N = 100$ and ${\Delta _v} = 42.05\;{{\rm m}^{- 1}}$. It is important to mention that $N = 100$ is the minimum, total number of ${\boldsymbol v}$-domain grid points (or pseudo-modes) because ${H_\alpha}$ is one-dimensional in ${\boldsymbol v}$, i.e., only a function of ${v_x}$.

With the ${\boldsymbol \rho}$ and ${\boldsymbol v}$ grids set, we are now ready to implement the simulation. The MATLAB script (.m file) that executes the simulation is provided in Code 1, Ref. [118]. The source code is explained line by line in the accompanying document (Ref. [118]).

 figure: Fig. 3.

Fig. 3. EGPSM source-plane CSDM results ${W_{\alpha \beta}}({{x_1},0,{x_2},0})$: (a)–(d) real and imaginary parts of ${W_{\textit{xx}}}$ theory and simulation, (e)–(h) real and imaginary parts of ${W_{\textit{xy}}}$ theory and simulation, (i)–(l) real and imaginary parts of ${W_{\textit{yx}}}$ theory and simulation, and (m)–(p) real and imaginary parts of ${W_{\textit{yy}}}$ theory and simulation.

Download Full Size | PDF

C. Results

1. Source-Plane Stokes Parameters

Figures 2–5 show the results of the simulation. In Fig. 2, we display the input- or source-plane Stokes parameters and speckle contrast. The first column shows the theory obtained from the EGPSM $\textbf{W}$ in Eq. (43)—note that the theoretical speckle contrast is

$${\cal C}({\boldsymbol \rho} ) = \sqrt {\frac{{1 + {{[{{\cal P}({\boldsymbol \rho} )} ]}^2}}}{2}} ,$$
where ${\cal P}$ is the degree of polarization [3,6,8,13]. The second column shows the simulated results after $T = 10,000$ trials (see [118]). Each row of Fig. 2 corresponds to a Stokes parameter, with the last row showing the speckle contrast. The theoretical and simulated results in each row are displayed using the same false color scale defined by the bar at row’s end. The rows and columns of the figure are labeled for the reader’s convenience.
 figure: Fig. 4.

Fig. 4. EGPSM source-plane field realizations: (a) and (b) $| {{E_x}} |$ and $\arg ({{E_x}})$, (c) and (d) $| {{E_y}} |$ and $\arg ({{E_y}})$, (e) and (f) $| {{E_x}/{\tau _x}} |$ and $\arg ({{E_x}/{\tau _x}})$, and (g) and (h) $| {{E_y}/{\tau _y}} |$ and $\arg ({{E_y}/{\tau _y}})$.

Download Full Size | PDF

 figure: Fig. 5.

Fig. 5. EGPSM observation-plane Stokes parameter results theory (solid lines) versus simulation ($\Delta$ symbols) for ${N_F} = \infty$, 15, 10, 5, and 1 (denoted by line color): (a) ${S_0}$, (b) ${S_1}$, (c) ${S_2}$, and (d) ${S_3}$.

Download Full Size | PDF

Overall, the agreement between theory and simulation is excellent. The maximum value of the simulated speckle contrast is approximately 1.02; therefore, the circular or ring-shaped artifacts evident in Fig. 2(j) are small in magnitude. The largest deviations from the theoretical ${\cal C}$ occur in regions of low light. Since ${\cal C} = {\sigma _I}/\langle I \rangle$ [8], errors in ${\sigma _I}$ are amplified in low-light regions resulting in “noisy” speckle contrast results. The reason the artifacts have circular symmetry is due to the ring-shaped EGPSM field realizations (discussed in more detail below) used to compute the simulated ${\cal C}$.

2. Source-Plane CSDM Results

The quality of the results in Fig. 2 proves that we are generating electric fields with the desired first-order (single-point, i.e., ${{\boldsymbol \rho}_1} = {{\boldsymbol \rho}_2} = {\boldsymbol \rho}$) behaviors. To validate that those electric fields possess the desired correlations (second-order field statistics), we must examine the CSDM results. Figure 3 displays images of the ${x_1}{-}{x_2}$ planar cuts through the elements of $\textbf{W}$. The layout of the figure is meant to resemble the $2 \times 2$ CSDM, where each “element” consists of four images arranged in a $2 \times 2$ pattern. The first column of each element displays the real (top row) and imaginary (bottom row) parts of ${W_{\alpha \beta}}$ theory [see Eq. (43)]; the second column shows the same for the simulated result. Both ${W_{\alpha \beta}}$ theory and simulation are encoded using the same color scale defined by the bars immediately to the right of the simulated images. The figure is labeled to aid the reader.

Again, we observe excellent agreement between the theoretical and simulated results. When considered together, Figs. 2 and 3 prove that we are indeed generating field realizations consistent with the desired EGPSM source. Note that the errors in the simulated imaginary parts of ${W_{\textit{xx}}}$ and ${W_{\textit{yy}}}$ [Figs. 3(d) and 3(p), respectively] are quite small (approximately ${10^{- 3}}$). In fact, they are well approximated by the speckle contrast (residual error) expression given in Eq. (42). To see this, let us consider the simulated ${W_{\textit{yy}}}$ results shown in Figs. 3(o) and 3(p). The maximum value of the real part of ${W_{\textit{yy}}}$ is approximately 0.1, while the error or “noise” in the imaginary part is roughly ${10^{- 3}}$. Thus, the residual or relative error is approximately 0.01. Recall that we computed the CSDM elements from $T = 10{,}000$ statistically independent ${\boldsymbol E}$. As a result, we should expect the residual error to be on the order of $1/\sqrt T \sim 0.01$, which is precisely what we observe. Performing the same analysis on the simulated ${W_{\textit{xx}}}$ results [Figs. 3(c) and 3(d)] yields a residual error on the order of 0.01—again consistent with the speckle contrast in Eq. (42).

3. EGPSM Field Realization

Figure 4 shows an example EGPSM field realization. Rows 1 and 2 display the magnitudes (column 1) and complex arguments (column 2) of ${E_x}$ and ${E_y}$, respectively. Rows 3 and 4 show the same, with the deterministic beam shapes [${\tau _x}$ and ${\tau _y}$; see Eq. (45)] effectively removed. Put simply, Figs. 4(e)–4(h) are the stochastic contributions to ${E_x}$ and ${E_y}$. Recall that ${E_x}$ and ${E_y}$ are thermal (CCG distributed) and, therefore, are fully developed speckle fields (i.e., ${| {{E_\alpha}} |^2}$ is exponentially distributed with a contrast equal to unity) [8,119–122]. They look nothing like classic speckle fields because their correlation functions [see Eq. (43)] are space (shift) variant. In fact, they both possess correlation functions that are space-invariant with respect to the radial coordinate $\rho$. Consequently, the speckles in Figs. 4(e)–4(h) are rotationally sheared, producing ring-shaped speckles.

4. Propagation Results

Lastly, Fig. 5 plots the theoretical and simulated (solid lines and $\Delta$ symbols, respectively) Stokes parameters after propagating ${N_F} = {\max}^2 ({{\sigma _\alpha}})k/({2z}) = \infty$, 15, 10, 5, and 1 (denoted by line color). The $\rho$ slices through ${S_0}$–${S_3}$ are along the $\phi$ at which ${S_0}$ theory obtains its maximum value (${\phi _{{\rm max}}}$ in the figure). For the five simulated propagation distances, ${\phi _{{\rm max}}} \approx {120^ \circ}$, 130°, 125°, 135°, and 135°, respectively. Note that ${N_F} = \infty$ is equivalent to $z = 0$; therefore, those traces (the blue lines) are the $\rho$ slices at ${\phi _{{\rm max}}} \approx {120^ \circ}$ through Figs. 2(a)–2(h).

Like in Figs. 2 and 3, the theoretical and simulated results are in excellent agreement. As mentioned above, one of the key features of EGPSM sources is that they “self-focus,” viz., their peak irradiance or spectral density occurs after propagating a near-field distance from the source. Figure 5(a) clearly shows this behavior, as the max ${S_0}$ for ${N_F} = 10$ and 5 are both greater than that for ${N_F} = \infty$.

4. CONCLUSION

In this tutorial, we reviewed the current techniques for simulating random electromagnetic fields consistent with a given CSDM $\textbf{W}$. We began with approaches based on the coherent-mode representation (or expansion) of $\textbf{W}$. We presented two ways in which coherent modes could be used to generate optical fields and discussed the pros and cons of coherent-mode-based simulation methods. Although very powerful (generating an optical field requires a relatively small number of modes), coherent modes are rarely used in wave optics simulations because so few are known.

We then discussed simulation approaches based on the superposition rule, which, in contrast to coherent modes, can be used to generate field realizations given any $\textbf{W}$. We presented two ways (analogous to coherent modes) in which the superposition rule could be applied or interpreted to synthesize optical fields. We explained how to properly sample and discretize field realizations by applying the Shannon–Nyquist sampling criterion and then discussed the strengths and weaknesses of using the superposition rule in random field simulations.

Lastly, we presented an example simulation where we generated (using the superposition rule) and propagated (through free space) an electromagnetic partially coherent beam from the literature. First, we performed the sampling analysis for the simulation, determining the maximum grid spacings and minimum numbers of grid points. Then, we analyzed the results by comparing the simulated statistics to theory. The agreement was excellent validating that we had indeed generated the desired random beam.

As statistical optics and optical coherence theory are applied in new technologies, simulation of stochastic electromagnetic fields becomes ever more important. The primary purpose of this tutorial is to aid scientists and engineers working on these new applications. Lastly, this tutorial is also meant to assist teachers of statistical optics. Using the techniques presented herein, instructors can create demonstrations of foundational statistical optics concepts, which augment the theory presented in the classic texts by Wolf, Mandel, and Goodman [6,13,20].

APPENDIX A: LIST OF ACRONYMS AND SYMBOLS

${{\boldsymbol \psi}_n}$Coherent-mode expansion vector eigenfunctions [see Eq. (7)]
${{\boldsymbol E}_t}$$2 \times 1$ transverse electric field vector
$\Delta$Spatial- or ${\boldsymbol \rho}$-domain sample spacing [see Eq. (36)]
${\Delta _v}$${\boldsymbol v}$-domain sample spacing [see Eq. (38)]
${\lambda _n}$Coherent-mode expansion eigenvalues [see Eq. (7)]
$\textbf{G}$Dyadic Green’s function or impulse response
$\textbf{H}$$2 \times 2$ diagonal kernel matrix [see Eq. (18)]
$\textbf{p}$$2 \times 2$ Hermitian, nonnegative definite matrix [see Eq. (18)]
$\textbf{W}$Cross-spectral density matrix
${\textbf{W}_t}$$2 \times 2$ transverse cross-spectral density matrix
${\cal C}$Speckle contrast
$\omega$Radian frequency
${\rm j}$$\sqrt {- 1}$
${D_p}$Maximum diameter of the ${\boldsymbol v}$-domain grid
${E_\alpha}$$\alpha = x,y$ Component of the electric field ${\boldsymbol E}$
${f_{{\rm max}}}$Maximum frequency of ${\boldsymbol E}$ in the spatial or ${\boldsymbol \rho}$ domain
${f_{v,{\rm max}}}$Maximum frequency of ${\boldsymbol E}$ in the ${\boldsymbol v}$ domain
$k$Wavenumber
$L$Physical side-length of the spatial- or ${\boldsymbol \rho}$-domain grid
$N$Minimum number of ${\boldsymbol v}$-domain grid points [see Eq. (38)]
$P$Minimum number of spatial- or ${\boldsymbol \rho}$-domain grid points [see Eq. (36)]
$T$Number of Monte Carlo trials [see Eq. (42)]
$W$Cross-spectral density function
${W_{\alpha \beta}}$$\alpha ,\beta = x,y$ Element of cross-spectral density matrix $\textbf{W}$
${H_\alpha}$$\alpha = x,y$ Component of diagonal kernel matrix $\textbf{H}$ [see Eq. (18)]
${p_{\alpha \beta}}$$\alpha ,\beta = x,y$ Element of $2 \times 2$ matrix $\textbf{p}$ [see Eq. (18)]
CCGCircular complex Gaussian
CSDCross-spectral density
CSDMCross-spectral density matrix
EGPSMElectromagnetic Gaussian pseudo-Schell model
EGSMElectromagnetic Gaussian Schell-model
FFTFast Fourier transform
GSMGAUSSIAN Schell-model
NUCNonuniformly correlated
WSSWide-sense stationary

Acknowledgment

The views expressed in this paper are those of the authors and do not reflect the official policy or position of the U.S. Air Force, the Department of Defense, or the U.S. government.

Disclosures

The author declares no conflicts of interest. This tutorial is a summary of Chapters 2 and 4 of the author’s forthcoming book Computational Optical Coherence and Statistical Optics.

Data availability

MATLAB files (.m files) related to the simulation are provided in Code 1, Ref. [118], along with a document that explains the source code.

REFERENCES

1. G. Gbur and T. Visser, “The structure of partially coherent fields,” in Progress in Optics, E. Wolf, ed. (Elsevier, 2010), Vol. 55, Chap. 5, pp. 285–341.

2. G. Gbur, “Partially coherent beam propagation in atmospheric turbulence,” J. Opt. Soc. Am. A 31, 2038–2045 (2014). [CrossRef]  

3. O. Korotkova, Random Light Beams: Theory and Applications (CRC Press, 2014).

4. Y. Cai, Y. Chen, J. Yu, X. Liu, and L. Liu, “Generation of partially coherent beams,” in Progress in Optics, T. D. Visser, ed. (Elsevier, 2017), Vol. 62, Chap. 3, pp. 157–223.

5. O. Korotkova and G. Gbur, “Applications of optical coherence theory,” in Progress in Optics, T. D. Visser, ed. (Elsevier, 2020), Vol. 65, Chap. 4, pp. 44–104.

6. J. W. Goodman, Statistical Optics, 2nd ed. (Wiley, 2015).

7. O. V. Angelsky, ed., Introduction to Singular Correlation Optics (SPIE, 2019).

8. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications, 2nd ed. (SPIE, 2020).

9. O. Korotkova, Theoretical Statistical Optics (World Scientific, 2021).

10. The Event Horizon Telescope Collaboration, “First M87 Event Horizon Telescope results. IV. Imaging the central supermassive black hole,” Astrophys. J. Lett. 875, L4 (2019). [CrossRef]  

11. The Event Horizon Telescope Collaboration, “First M87 Event Horizon Telescope results. VI. The shadow and mass of the central black hole,” Astrophys. J. Lett. 875, L6 (2019). [CrossRef]  

12. E. Wolf, “Unified theory of coherence and polarization of random electromagnetic beams,” Phys. Lett. A 312, 263–267 (2003). [CrossRef]  

13. E. Wolf, Introduction to the Theory of Coherence and Polarization of Light (Cambridge University, 2007).

14. T. Voipio, T. Setälä, and A. T. Friberg, “Partial polarization theory of pulsed optical beams,” J. Opt. Soc. Am. A 30, 71–81 (2013). [CrossRef]  

15. T. Voipio, T. Setälä, and A. T. Friberg, “Coherent-mode representation of partially polarized pulsed electromagnetic beams,” J. Opt. Soc. Am. A 30, 2433–2443 (2013). [CrossRef]  

16. A. T. Friberg and T. Setälä, “Electromagnetic theory of optical coherence,” J. Opt. Soc. Am. A 33, 2431–2442 (2016). [CrossRef]  

17. S. A. Collins, “Lens-system diffraction integral written in terms of matrix optics,” J. Opt. Soc. Am. 60, 1168–1177 (1970). [CrossRef]  

18. J. W. Goodman, Introduction to Fourier Optics, 4th ed. (W. H. Freeman and Company, 2017).

19. J. D. Gaskill, Linear Systems, Fourier Transforms, and Optics (Wiley, 1978).

20. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University, 1995).

21. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, 7th ed. (Cambridge University, 1999).

22. T. Hansen and A. Yaghjian, Plane-Wave Theory of Time-Domain Fields (Near-Field Scanning Applications) (IEEE, 1999).

23. E. J. Rothwell and M. J. Cloud, Electromagnetics, 2nd ed. (CRC Press, 2008).

24. A. Yariv and P. Yeh, Optical Waves in Crystals: Propagation and Control of Laser Radiation (Wiley, 2003).

25. T. G. Mackay and A. Lakhtakia, Electromagnetic Anisotropy and Bianisotropy: A Field Guide (World Scientific, 2010).

26. W. Chew, Waves and Fields in Inhomogeneous Media (IEEE, 1999).

27. A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE, 1999).

28. L. C. Andrews and R. L. Phillips, Laser Beam Propagation through Random Media, 2nd ed. (SPIE, 2005).

29. T. Magalhães and J. M. Rebordão, “Numerical simulations of spectral shifts in the far-field spectrum of light due to source correlations,” J. Opt. Soc. Am. A 35, 423–430 (2018). [CrossRef]  

30. T. C. Magalhães and J. M. Rebordão, “Simulation of partially coherent light propagation using parallel computing devices,” Proc. SPIE 10453, 439–451 (2017). [CrossRef]  

31. T. E. C. Magalhães and J. M. Rebordão, “PyWolf: a PyOpenCL implementation for simulating the propagation of partially coherent light,” Comput. Phys. Commun. 276, 108336 (2022). [CrossRef]  

32. F. Gori, “Mode propagation of the field generated by Collett-Wolf Schell-model sources,” Opt. Commun. 46, 149–154 (1983). [CrossRef]  

33. J. Huttunen, A. T. Friberg, and J. Turunen, “Scattering of partially coherent electromagnetic fields by microstructured media,” Phys. Rev. E 52, 3081–3092 (1995). [CrossRef]  

34. T. Shirai, A. Dogariu, and E. Wolf, “Mode analysis of spreading of partially coherent beams propagating through atmospheric turbulence,” J. Opt. Soc. Am. A 20, 1094–1102 (2003). [CrossRef]  

35. J. Lindberg, T. Setälä, M. Kaivola, and A. T. Friberg, “Spatial coherence effects in light scattering from metallic nanocylinders,” J. Opt. Soc. Am. A 23, 1349–1358 (2006). [CrossRef]  

36. A. S. Ostrovsky, Coherent-Mode Representations in Optics (SPIE, 2006).

37. P. Vahimaa and J. Turunen, “Finite-elementary-source model for partially coherent radiation,” Opt. Express 14, 1376–1381 (2006). [CrossRef]  

38. A. Burvall, A. Smith, and C. Dainty, “Elementary functions: propagation of partially coherent light,” J. Opt. Soc. Am. A 26, 1721–1729 (2009). [CrossRef]  

39. Y. Gu and G. Gbur, “Scintillation of pseudo-Bessel correlated beams in atmospheric turbulence,” J. Opt. Soc. Am. A 27, 2621–2629 (2010). [CrossRef]  

40. J. Turunen, “Elementary-field representations in partially coherent optics,” J. Mod. Opt. 58, 509–527 (2011). [CrossRef]  

41. M. Singh, J. Tervo, and J. Turunen, “Elementary-field analysis of partially coherent beam shaping,” J. Opt. Soc. Am. A 30, 2611–2617 (2013). [CrossRef]  

42. A. Smith and C. Dainty, “Numerical modeling of spatial coherence using the elementary function method,” Appl. Opt. 52, 2815–2827 (2013). [CrossRef]  

43. Y. Gu and G. Gbur, “Scintillation of nonuniformly correlated beams in atmospheric turbulence,” Opt. Lett. 38, 1395–1397 (2013). [CrossRef]  

44. M. Singh, H. Lajunen, J. Tervo, and J. Turunen, “Imaging with partially coherent light: elementary-field approach,” Opt. Express 23, 28132–28140 (2015). [CrossRef]  

45. P. Ma, B. Kacerovská, R. Khosravi, C. Liang, J. Zeng, X. Peng, C. Mi, Y. E. Monfared, Y. Zhang, F. Wang, and Y. Cai, “Numerical approach for studying the evolution of the degrees of coherence of partially coherent beams propagation through an ABCD optical system,” Appl. Sci. 9, 2084 (2019). [CrossRef]  

46. M. W. Hyde, “Direct synthesis of partially coherent fields at the outputs of ABCD optical systems,” in IEEE Aerospace Conference (2021), pp. 1–6.

47. M. Aviñoá, R. Martínez-Herrero, and A. Carnicer, “Efficient calculation of highly focused electromagnetic Schell-model beams,” Opt. Express 29, 26220–26232 (2021). [CrossRef]  

48. D. G. Voelz, Computational Fourier Optics: A MATLAB Tutorial (SPIE, 2011).

49. J. D. Schmidt, Numerical Simulation of Optical Wave Propagation with Examples in MATLAB (SPIE, 2010).

50. S. W. Teare, Optics Using MATLAB (SPIE, 2017).

51. J. M. Jarem and P. P. Banerjee, Computational Methods for Electromagnetic and Optical Systems, 2nd ed. (CRC Press, 2011).

52. M. S. Wartak, Computational Photonics: An Introduction with MATLAB (Cambridge University, 2013).

53. P. De Santis, F. Gori, G. Guattari, and C. Palma, “Synthesis of partially coherent fields,” J. Opt. Soc. Am. A 3, 1258–1262 (1986). [CrossRef]  

54. T. Shirai, O. Korotkova, and E. Wolf, “A method of generating electromagnetic Gaussian Schell-model beams,” J. Opt. A 7, 232–237 (2005). [CrossRef]  

55. G. Piquero, F. Gori, P. Romanini, M. Santarsiero, R. Borghi, and A. Mondello, “Synthesis of partially polarized Gaussian Schell-model sources,” Opt. Commun. 208, 9–16 (2002). [CrossRef]  

56. G. J. Gbur, “Simulating fields of arbitrary spatial and temporal coherence,” Opt. Express 14, 7567–7578 (2006). [CrossRef]  

57. C. Rydberg and J. Bengtsson, “Efficient numerical representation of the optical field for the propagation of partially coherent radiation with a specified spatial and temporal coherence function,” J. Opt. Soc. Am. A 23, 1616–1625 (2006). [CrossRef]  

58. B. J. Davis, “Simulation of vector fields with arbitrary second-order correlations,” Opt. Express 15, 2837–2846 (2007). [CrossRef]  

59. A. S. Ostrovsky, G. Rodríguez-Zurita, C. Meneses-Fabián, M. A. Olvera-Santamaría, and C. Rickenstorff-Parrao, “Experimental generating the partially coherent and partially polarized electromagnetic source,” Opt. Express 18, 12864–12871 (2010). [CrossRef]  

60. D. Voelz, X. Xiao, and O. Korotkova, “Numerical modeling of Schell-model beams with arbitrary far-field patterns,” Opt. Lett. 40, 352–355 (2015). [CrossRef]  

61. M. Santarsiero, R. Martínez-Herrero, D. Maluenda, J. C. G. de Sande, G. Piquero, and F. Gori, “Synthesis of circularly coherent sources,” Opt. Lett. 42, 4115–4118 (2017). [CrossRef]  

62. Y. Cai, Y. Chen, and F. Wang, “Generation and propagation of partially coherent beams with nonconventional correlation functions: a review,” J. Opt. Soc. Am. A 31, 2083–2096 (2014). [CrossRef]  

63. X. Chen, J. Li, S. M. H. Rafsanjani, and O. Korotkova, “Synthesis of Im-Bessel correlated beams via coherent modes,” Opt. Lett. 43, 3590–3593 (2018). [CrossRef]  

64. F. Wang, H. Lv, Y. Chen, Y. Cai, and O. Korotkova, “Three modal decompositions of Gaussian Schell-model sources: comparative analysis,” Opt. Express 29, 29676–29689 (2021). [CrossRef]  

65. M. W. Hyde, “Stochastic complex transmittance screens for synthesizing general partially coherent sources,” J. Opt. Soc. Am. A 37, 257–264 (2020). [CrossRef]  

66. M. W. Hyde IV, “Synthesizing general electromagnetic partially coherent sources from random, correlated complex screens,” Optics 1, 97–113 (2020). [CrossRef]  

67. M. W. Hyde and O. Korotkova, “Pseudo-modal expansions for generating random electromagnetic beams,” J. Opt. Soc. Am. A 39, 545–551 (2022). [CrossRef]  

68. J. Tervo, T. Setälä, and A. T. Friberg, “Degree of coherence for electromagnetic fields,” Opt. Express 11, 1137–1143 (2003). [CrossRef]  

69. T. Setälä, J. Tervo, and A. T. Friberg, “Theorems on complete electromagnetic coherence in the space-time domain,” Opt. Commun. 238, 229–236 (2004). [CrossRef]  

70. E. Wolf, “New spectral representation of random sources and of the partially coherent fields that they generate,” Opt. Commun. 38, 3–6 (1981). [CrossRef]  

71. E. Wolf, “New theory of partial coherence in the space-frequency domain. Part I: spectra and cross spectra of steady-state sources,” J. Opt. Soc. Am. 72, 343–351 (1982). [CrossRef]  

72. F. Gori, M. Santarsiero, R. Simon, G. Piquero, R. Borghi, and G. Guattari, “Coherent-mode decomposition of partially polarized, partially coherent sources,” J. Opt. Soc. Am. A 20, 78–84 (2003). [CrossRef]  

73. J. Tervo, T. Setälä, and A. T. Friberg, “Theory of partially coherent electromagnetic fields in the space-frequency domain,” J. Opt. Soc. Am. A 21, 2205–2215 (2004). [CrossRef]  

74. H. Lajunen, J. Tervo, and P. Vahimaa, “Overall coherence and coherent-mode expansion of spectrally partially coherent plane-wave pulses,” J. Opt. Soc. Am. A 21, 2117–2123 (2004). [CrossRef]  

75. H. Lajunen, P. Vahimaa, and J. Tervo, “Theory of spatially and spectrally partially coherent pulses,” J. Opt. Soc. Am. A 22, 1536–1545 (2005). [CrossRef]  

76. T. Voipio, K. Blomstedt, T. Setälä, and A. T. Friberg, “Conservation of electromagnetic coherent-mode structure on propagation,” Opt. Commun. 340, 93–101 (2015). [CrossRef]  

77. B. Rodenburg, M. Mirhosseini, O. S. Magaña-Loaiza, and R. W. Boyd, “Experimental generation of an optical field with arbitrary spatial coherence properties,” J. Opt. Soc. Am. B 31, A51–A55 (2014). [CrossRef]  

78. X. Zhu, F. Wang, C. Zhao, Y. Cai, and S. A. Ponomarenko, “Experimental realization of dark and antidark diffraction-free beams,” Opt. Lett. 44, 2260–2263 (2019). [CrossRef]  

79. A. Bhattacharjee, R. Sahu, and A. K. Jha, “Generation of a Gaussian Schell-model field as a mixture of its coherent modes,” J. Opt. 21, 105601 (2019). [CrossRef]  

80. F. Gori, M. Santarsiero, and R. Borghi, “Modal expansion forJ0-correlated electromagnetic sources,” Opt. Lett. 33, 1857–1859 (2008). [CrossRef]  

81. K. Kim and E. Wolf, “A scalar-mode representation of stochastic, planar, electromagnetic sources,” Opt. Commun. 261, 19–22 (2006). [CrossRef]  

82. K. Kim, “Vector and scalar modes in coherent mode representation of electromagnetic beams,” J. Opt. Soc. Korea 12, 103–106 (2008). [CrossRef]  

83. D. F. V. James, “Change of polarization of light beams on propagation in free space,” J. Opt. Soc. Am. A 11, 1641–1643 (1994). [CrossRef]  

84. O. Korotkova, M. Salem, and E. Wolf, “Beam conditions for radiation generated by an electromagnetic Gaussian Schell-model source,” Opt. Lett. 29, 1173–1175 (2004). [CrossRef]  

85. O. Korotkova and E. Wolf, “Changes in the state of polarization of a random electromagnetic beam on propagation,” Opt. Commun. 246, 35–43 (2005). [CrossRef]  

86. F. Gori, M. Santarsiero, R. Borghi, and V. Ramírez-Sánchez, “Realizability condition for electromagnetic Schell-model sources,” J. Opt. Soc. Am. A 25, 1016–1021 (2008). [CrossRef]  

87. H. Roychowdhury and O. Korotkova, “Realizability conditions for electromagnetic Gaussian Schell-model sources,” Opt. Commun. 249, 379–385 (2005). [CrossRef]  

88. B. J. Davis and R. W. Schoonover, “Computationally efficient coherent-mode representations,” Opt. Lett. 34, 923–925 (2009). [CrossRef]  

89. A. Schell, “A technique for the determination of the radiation pattern of a partially coherent aperture,” IEEE Trans. Antennas Propag. 15, 187–188 (1967). [CrossRef]  

90. F. Gori and M. Santarsiero, “Devising genuine spatial correlation functions,” Opt. Lett. 32, 3531–3533 (2007). [CrossRef]  

91. F. Gori, V. Ramírez-Sánchez, M. Santarsiero, and T. Shirai, “On genuine cross-spectral density matrices,” J. Opt. A 11, 085706 (2009). [CrossRef]  

92. T. Shirai, “Interpretation of the recipe for synthesizing genuine cross-spectral density matrices,” Opt. Commun. 283, 4478–4483 (2010). [CrossRef]  

93. Y. Chen, F. Wang, and Y. Cai, “Partially coherent light beam shaping via complex spatial coherence structure engineering,” Adv. Phys. X 7, 2009742 (2022). [CrossRef]  

94. R. Martínez-Herrero, P. M. Mejías, and F. Gori, “Genuine cross-spectral densities and pseudo-modal expansions,” Opt. Lett. 34, 1399–1401 (2009). [CrossRef]  

95. R. Martínez-Herrero and P. M. Mejías, “Elementary-field expansions of genuine cross-spectral density matrices,” Opt. Lett. 34, 2303–2305 (2009). [CrossRef]  

96. J. Tervo, J. Turunen, P. Vahimaa, and F. Wyrowski, “Shifted-elementary-mode representation for partially coherent vectorial fields,” J. Opt. Soc. Am. A 27, 2004–2014 (2010). [CrossRef]  

97. E. Wolf, “Coherence properties of partially polarized electromagnetic radiation,” Nuovo Cimento 13, 1165–1181 (1959). [CrossRef]  

98. E. L. O’Neill, Introduction to Statistical Optics (Dover, 1992).

99. C. Brosseau, Fundamentals of Polarized Light: A Statistical Optics Approach (Wiley, 1998).

100. S. M. Kay, Intuitive Probability and Random Processes Using MATLAB (Springer, 2006).

101. D. S. Watkins, Fundamentals of Matrix Computations, 2nd ed. (Wiley, 2002).

102. A. Jerri, “The Shannon sampling theorem—its various extensions and applications: a tutorial review,” Proc. IEEE 65, 1565–1596 (1977). [CrossRef]  

103. R. J. Marks, J. F. Walkup, and M. O. Hagler, “A sampling theorem for space-variant systems,” J. Opt. Soc. Am. 66, 918–921 (1976). [CrossRef]  

104. R. J. Marks, “Sampling theory for linear integral transforms,” Opt. Lett. 6, 7–9 (1981). [CrossRef]  

105. R. Marks, J. Walkup, and M. Hagler, “Sampling theorems for linear shift-variant systems,” IEEE Trans. Circuits Syst. 25, 228–233 (1978). [CrossRef]  

106. D. G. Voelz and M. C. Roggemann, “Digital simulation of scalar optical diffraction: revisiting chirp function sampling criteria and consequences,” Appl. Opt. 48, 6132–6142 (2009). [CrossRef]  

107. D. P. Kelly, “Numerical calculation of the Fresnel transform,” J. Opt. Soc. Am. A 31, 755–764 (2014). [CrossRef]  

108. W. Zhang, H. Zhang, and G. Jin, “Frequency sampling strategy for numerical diffraction calculations,” Opt. Express 28, 39916–39932 (2020). [CrossRef]  

109. F. Gori, “Collett-Wolf sources and multimode lasers,” Opt. Commun. 34, 301–305 (1980). [CrossRef]  

110. A. Starikov and E. Wolf, “Coherent-mode representation of Gaussian Schell-model sources and of their radiation fields,” J. Opt. Soc. Am. 72, 923–928 (1982). [CrossRef]  

111. A. Starikov, “Effective number of degrees of freedom of partially coherent sources,” J. Opt. Soc. Am. 72, 1538–1544 (1982). [CrossRef]  

112. J. C. G. de Sande, R. Martínez-Herrero, G. Piquero, M. Santarsiero, and F. Gori, “Pseudo-Schell model sources,” Opt. Express 27, 3963–3977 (2019). [CrossRef]  

113. R. Martínez-Herrero, D. Maluenda, G. Piquero, J. C. G. de Sande, M. Santarsiero, and F. Gori, “Partially polarized pseudo-Schell model sources,” Proc. SPIE 10453, 398–403 (2017). [CrossRef]  

114. R. Martínez-Herrero, G. Piquero, J. C. González de Sande, M. Santarsiero, and F. Gori, “Besinc pseudo-Schell model sources with circular coherence,” Appl. Sci. 9, 2716 (2019). [CrossRef]  

115. M. Santarsiero, R. Martínez-Herrero, G. Piquero, J. C. G. de Sande, and F. Gori, “Modal analysis of pseudo-Schell model sources,” Photonics 8, 449 (2021). [CrossRef]  

116. R. Martínez-Herrero, G. Piquero, J. C. González de Sande, M. Santarsiero, and F. Gori, “Propagation features of a class of beams radiated from pseudo-Schell model vectorial sources,” EPJ Web Conf. 255, 12012 (2021). [CrossRef]  

117. R. Martínez-Herrero, G. Piquero, M. Santarsiero, F. Gori, and J. C. González de Sande, “A class of vectorial pseudo-Schell model sources with structured coherence and polarization,” Opt. Laser Technol. 152, 108079 (2022). [CrossRef]  

118. M. Hyde, “MATLAB R2018b electromagnetic Gaussian pseudo-Schell model beam simulation,” figshare, 2022, https://doi.org/10.6084/m9.figshare.19918237 [retrieved Sept. 27, 2022].

119. J. Goodman, “Statistical properties of laser speckle patterns,” in Laser Speckle and Related Phenomena, J. C. Dainty, ed., Vol. 9 of Topics in Applied Physics (Springer-Verlag, 1975), pp. 9–75.

120. J. Dainty, “The statistics of speckle patterns,” in Progress in Optics, E. Wolf, ed. (Elsevier, 1977), Vol. 14, Chap. 1, pp. 1–46.

121. I. Freund, “‘1001’ correlations in random wave fields,” Waves Random Media 8, 119–158 (1998). [CrossRef]  

122. G. J. Gbur, Singular Optics (CRC Press, 2016).

Supplementary Material (1)

NameDescription
Code 1       MATLAB scripts (.m files) associated with "Simulating random optics fields: tutorial"

Data availability

MATLAB files (.m files) related to the simulation are provided in Code 1, Ref. [118], along with a document that explains the source code.

118. M. Hyde, “MATLAB R2018b electromagnetic Gaussian pseudo-Schell model beam simulation,” figshare, 2022, https://doi.org/10.6084/m9.figshare.19918237 [retrieved Sept. 27, 2022].

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (5)
Fig. 1.
Fig. 1. Summary of approaches for simulating random optical fields. The green and red arrows show the directions of the pros (${+}$) and cons (${-}$) comparisons. For instance, the top arrow points to the coherent-mode and superposition-rule boxes. Thus, the pros and cons listed in those boxes relate coherent-mode-based and superposition-rule-based simulation approaches. Likewise, the bottom arrow points to methods 1 and 2 in both the coherent-mode and superposition-rule boxes. Therefore, the pros and cons listed under those methods relate coherent-mode method 1 to coherent-mode method 2 and superposition-rule method 1 to superposition-rule method 2.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. EGPSM source-plane Stokes parameter and speckle contrast ${\cal C}$ results: (a) and (b) ${S_0}$ theory and simulation, (c) and (d) ${S_1}$ theory and simulation, (e) and (f) ${S_2}$ theory and simulation, (g) and (h) ${S_3}$ theory and simulation, and (i) and (j) ${\cal C}$ theory and simulation.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. EGPSM source-plane CSDM results ${W_{\alpha \beta}}({{x_1},0,{x_2},0})$: (a)–(d) real and imaginary parts of ${W_{\textit{xx}}}$ theory and simulation, (e)–(h) real and imaginary parts of ${W_{\textit{xy}}}$ theory and simulation, (i)–(l) real and imaginary parts of ${W_{\textit{yx}}}$ theory and simulation, and (m)–(p) real and imaginary parts of ${W_{\textit{yy}}}$ theory and simulation.

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. EGPSM source-plane field realizations: (a) and (b) $| {{E_x}} |$ and $\arg ({{E_x}})$, (c) and (d) $| {{E_y}} |$ and $\arg ({{E_y}})$, (e) and (f) $| {{E_x}/{\tau _x}} |$ and $\arg ({{E_x}/{\tau _x}})$, and (g) and (h) $| {{E_y}/{\tau _y}} |$ and $\arg ({{E_y}/{\tau _y}})$.

View in Article | Download Full Size | PDF

Fig. 5.
Fig. 5. EGPSM observation-plane Stokes parameter results theory (solid lines) versus simulation ($\Delta$ symbols) for ${N_F} = \infty$, 15, 10, 5, and 1 (denoted by line color): (a) ${S_0}$, (b) ${S_1}$, (c) ${S_2}$, and (d) ${S_3}$.

View in Article | Download Full Size | PDF

Equations (49)

Equations on this page are rendered with MathJax. Learn more.

W ( r 1 , r 2 , ω ) = ⨌ − ∞ ∞ G ( r 1 , ω ; ρ 1 ′ ) W t ( ρ 1 ′ , ρ 2 ′ , ω ) × G † ( r 2 , ω ; ρ 2 ′ ) d 2 ρ 1 ′ d 2 ρ 2 ′ ,
W t ( ρ 1 ′ , ρ 2 ′ , ω ) = ⟨ [ E x ( ρ 1 ′ , ω ) E y ( ρ 1 ′ , ω ) ] [ E x ( ρ 2 ′ , ω ) E y ( ρ 2 ′ , ω ) ] † ⟩ .
E ( r , ω ) = ∬ − ∞ ∞ G ( r , ω ; ρ ′ ) E t ( ρ ′ , ω ) d 2 ρ ′ ,
∫ − ∞ ∞ ⨌ − ∞ ∞ | W α β ( ρ 1 , ρ 2 , ω ) | 2 d 2 ρ 1 d 2 ρ 2 d ω < ∞ ,
W ( ρ 1 , ρ 2 , ω ) = W † ( ρ 2 , ρ 1 , ω ) ,
⨌ − ∞ ∞ g † ( ρ 2 , ω ) W ( ρ 1 , ρ 2 , ω ) g ( ρ 1 , ω ) d 2 ρ 1 d 2 ρ 2 ≥ 0 ,
W ( ρ 1 , ρ 2 ) = ∑ n λ n ψ n ( ρ 1 ) ψ n † ( ρ 2 ) ,
∬ − ∞ ∞ ψ m † ( ρ ) ψ n ( ρ ) d 2 ρ = δ mn .
∬ − ∞ ∞ W ( ρ 1 , ρ 2 ) ψ n ( ρ 2 ) d 2 ρ 2 = λ n ψ n ( ρ 1 ) .
E n ( ρ ) = λ n ψ n ( ρ ) ,
E ( ρ ) = ∑ n a n ψ n ( ρ ) ,
⟨ E ( ρ 1 ) E † ( ρ 2 ) ⟩ = ∑ n ∑ m ⟨ a n a m ∗ ⟩ ψ n ( ρ 1 ) ψ m † ( ρ 2 ) ,
⟨ ( a n r ) 2 ⟩ = ⟨ ( a n i ) 2 ⟩ = λ n / 2.
E ( ρ ) = x ^ ∑ n a n ψ n x ( ρ ) + y ^ ∑ m b m ψ m y ( ρ ) ,
⟨ a n r ⟩ = ⟨ a n i ⟩ = ⟨ b m r ⟩ = ⟨ b m i ⟩ = 0 , ⟨ ( a n r ) 2 ⟩ = ⟨ ( a n i ) 2 ⟩ = λ n x / 2 , ⟨ ( b m r ) 2 ⟩ = ⟨ ( b m i ) 2 ⟩ = λ m y / 2 ,
⟨ a n r b m r ⟩ = ⟨ a n i b m i ⟩ = R nm r / 2 , ⟨ a n i b m r ⟩ = − ⟨ a n r b m i ⟩ = R nm i / 2 ,
R nm = ⨌ − ∞ ∞ W xy ( ρ 1 , ρ 2 ) [ ψ n x ( ρ 1 ) ] ∗ ψ m y ( ρ 2 ) d 2 ρ 1 d 2 ρ 2 .
W ( ρ 1 , ρ 2 ) = ∬ − ∞ ∞ H ( ρ 1 ; v ) p ( v ) H † ( ρ 2 ; v , ) d 2 v ,
p ( v ) = [ p xx ( v ) p xy ( v ) p xy ∗ ( v ) p yy ( v ) ] ,
H ( ρ ; v ) = [ H x ( ρ ; v ) 0 0 H y ( ρ ; v ) ] ,
p ( v ) = [ A ( v ) 0 0 A ( v ) ] + [ B ( v ) C ( v ) C ∗ ( v ) D ( v ) ] = [ A ( v ) 0 0 A ( v ) ] + [ u x ( v ) u y ( v ) ] [ u x ( v ) u y ( v ) ] † ,
u x ( v ) = B ( v ) exp ⁡ { j 1 2 arg ⁡ [ C ( v ) ] } , u y ( v ) = D ( v ) exp ⁡ { − j 1 2 arg ⁡ [ C ( v ) ] } .
W ( ρ 1 , ρ 2 ) = ∬ − ∞ ∞ A ( v ) [ H x ( ρ 1 ; v ) 0 ] [ H x ( ρ 2 ; v ) 0 ] † d 2 v + ∬ − ∞ ∞ A ( v ) [ 0 H y ( ρ 1 ; v ) ] [ 0 H y ( ρ 2 ; v ) ] † d 2 v + ∬ − ∞ ∞ [ u x ( v ) H x ( ρ 1 ; v ) u y ( v ) H y ( ρ 1 ; v ) ] [ u x ( v ) H x ( ρ 2 ; v ) u y ( v ) H y ( ρ 2 ; v ) ] † d 2 v .
E 1 ( ρ ; v ) = x ^ A ( v ) H x ( ρ ; v ) , E 2 ( ρ ; v ) = y ^ A ( v ) H y ( ρ ; v ) , E 3 ( ρ ; v ) = x ^ u x ( v ) H x ( ρ ; v ) + y ^ u y ( v ) H y ( ρ ; v ) = x ^ B ( v ) exp ⁡ { j 1 2 arg ⁡ [ C ( v ) ] } H x ( ρ ; v ) + y ^ D ( v ) exp { − j 1 2 arg ⁡ [ C ( v ) ] } H y ( ρ ; v ) .
A ( v ) = 1 2 { [ p xx ( v ) + p yy ( v ) ] − [ p xx ( v ) − p yy ( v ) ] 2 + 4 | p xy ( v ) | 2 } , B ( v ) = 1 2 { [ p xx ( v ) − p yy ( v ) ] + [ p xx ( v ) − p yy ( v ) ] 2 + 4 | p xy ( v ) | 2 } , C ( v ) = p xy ( v ) , D ( v ) = 1 2 { [ p yy ( v ) − p xx ( v ) ] + [ p xx ( v ) − p yy ( v ) ] 2 + 4 | p xy ( v ) | 2 } .
W α β ( ρ 1 , ρ 2 ) = ∬ − ∞ ∞ p α β ( v ) H α ( ρ 1 ; v ) H β ∗ ( ρ 2 ; v ) d 2 v .
E α ( ρ ) = ∬ − ∞ ∞ r α ( v ) [ 1 2 p α α ( v ) ] 1 / 2 H α ( ρ ; v ) d 2 v ,
⟨ E x ( ρ 1 ) E y ∗ ( ρ 2 ) ⟩ = ⨌ − ∞ ∞ ⟨ r x ( v 1 ) r y ∗ ( v 2 ) ⟩ × 1 2 [ p xx ( v 1 ) p yy ( v 2 ) ] 1 / 2 H x ( ρ 1 ; v 1 ) H y ∗ ( ρ 2 ; v 2 ) d 2 v 1 d 2 v 2 .
⟨ r x ( v 1 ) r y ∗ ( v 2 ) ⟩ = 2 p xy ( v 1 ) [ p xx ( v 1 ) p yy ( v 2 ) ] 1 / 2 δ ( v 1 − v 2 ) ,
R r ( v ) = ⟨ r x r ( v ) r y r ( v ) ⟩ = ⟨ r x i ( v ) r y i ( v ) ⟩ = R e [ p xy ( v ) ] [ p xx ( v ) p yy ( v ) ] 1 / 2 , R i ( v ) = ⟨ r x i ( v ) r y r ( v ) ⟩ = − ⟨ r x r v r y i ( v ) ⟩ = I m [ p xy ( v ) ] [ p xx ( v ) p yy ( v ) ] 1 / 2 .
E x ( ρ ) = ∬ − ∞ ∞ r x ( v ) [ 1 2 p xx ( v ) ] 1 / 2 H x ( ρ ; v ) d 2 v , E y ( ρ ) = ∬ − ∞ ∞ r y ( v ) [ 1 2 p yy ( v ) ] 1 / 2 H y ( ρ ; v ) d 2 v ,
Σ = [ ⟨ ( r x r ) 2 ⟩ ⟨ r x r r x i ⟩ ⟨ r x r r y r ⟩ ⟨ r x r r y i ⟩ ⟨ r x i r x r ⟩ ⟨ ( r x i ) 2 ⟩ ⟨ r x i r y r ⟩ ⟨ r x i r y i ⟩ ⟨ r y r r x r ⟩ ⟨ r y r r x i ⟩ ⟨ ( r y r ) 2 ⟩ ⟨ r y r r y i ⟩ ⟨ r y i r x r ⟩ ⟨ r y i r x i ⟩ ⟨ r y i r y r ⟩ ⟨ ( r y i ) 2 ⟩ ] = [ 1 0 R r − R i 0 1 R i R r R r R i 1 0 − R i R r 0 1 ] .
r x ( v ) = r 1 ( v ) + j r 2 ( v ) , r y ( v ) = r x ( v ) [ R r ( v ) − j R i ( v ) ] + [ r 3 ( v ) + j r 4 ( v ) ] × 1 − [ R r ( v ) ] 2 − [ R i ( v ) ] 2 ,
Δ ≤ 1 2 f m a x ,
H ~ α ( f ; v ) = ∬ − ∞ ∞ H α ( ρ ; v ) exp ⁡ ( − j 2 π f ⋅ ρ ) d 2 ρ .
P ≥ c e i l ( L / Δ ) = c e i l ( 2 L f m a x ) ,
f v , m a x = f v , m a x p + f v , m a x H ,
N ≥ c e i l ( D p / Δ v ) = c e i l ( 2 D p f v , m a x ) ,
Δ v = 2 π L = 2 π P Δ = 2 π N Δ , Δ = 2 π L v = 2 π N Δ v = 2 π P Δ v ,
∬ − ∞ ∞ G ( r ; ρ ′ ) H ( ρ ′ ; v ) d 2 ρ ′ ,
f m a x = f m a x G + f m a x H ,
C = 1 T .
W α β ( ρ 1 , ρ 2 ) = A α ρ 1 σ α cos ⁡ ( ϕ 1 − θ α ) exp ⁡ ( − ρ 1 2 σ α 2 ) × A β ρ 2 σ β cos ⁡ ( ϕ 2 − θ β ) exp ⁡ ( − ρ 2 2 σ β 2 ) B α β × exp ⁡ [ − ( ρ 1 − ρ 2 ) 2 δ α β 2 ] ,
B xx = B yy = 1 , B xy = B yx ∗ , | B xy | ≤ 1 , δ xy = δ yx , δ xx 2 + δ yy 2 2 ≤ δ xy ≤ δ xx δ yy | B xy | .
p α β ( v ) = B α β δ α β 2 4 π exp ⁡ ( − δ α β 2 4 v 2 ) , H α ( ρ ; v ) = τ α ( ρ ) h ( ρ ; v ) = A α ρ σ α cos ⁡ ( ϕ − θ α ) exp ( − ρ 2 σ α 2 ) exp ⁡ ( j ρ v x ) .
E α ( ρ , z ) = exp ⁡ ( j k z ) j λ z exp ⁡ ( j k 2 z ρ 2 ) ∬ − ∞ ∞ E α ( ρ ′ , 0 ) × exp ⁡ ( j k 2 z ρ ′ 2 ) exp ⁡ ( − j k z ρ ⋅ ρ ′ ) d 2 ρ ′ ,
Δ ≤ λ z D i n + 2 λ z f m a x , N ≥ D i n + D o u t Δ .
Δ v ≤ 2 π 2 max ( δ α β ) − ln ⁡ ϵ + D i n , N ≥ c e i l { 4 π ( − ln ⁡ ϵ ) max ( δ α β ) min ( δ α β ) [ 1 + D i n 2 max ( δ α β ) − ln ⁡ ϵ ] } .
C ( ρ ) = 1 + [ P ( ρ ) ] 2 2 ,
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use