Controlling the Spatial Coherence of an Optical Source Using a Spatial Filter

Featured Application

This technique for controlling spatial coherence can be used for beam shaping and reducing scintillation (or fades) in directed energy and free-space optical communication applications.

Abstract

This paper presents the theory for controlling the spectral degree of coherence via spatial filtering. Starting with a quasi-homogeneous partially coherent source, the cross-spectral density function of the field at the output of the spatial filter is found by applying Fourier and statistical optics theory. The key relation obtained from this analysis is a closed-form expression for the filter function in terms of the desired output spectral degree of coherence. This theory is verified with Monte Carlo wave-optics simulations of spatial coherence control and beam shaping for potential use in free-space optical communications and directed energy applications. The simulated results are found to be in good agreement with the developed theory. The technique presented in this paper will be useful in applications where coherence control is advantageous, e.g., directed energy, free-space optical communications, remote sensing, medicine, and manufacturing.

1. Introduction

Following the nomenclature and notation of Emil Wolf [1], a partially coherent source (PCS) is defined by its autocorrelation function, better known as the mutual coherence function (MCF):

$$\Gamma \left({\mathit{\rho }}_1,{\mathit{\rho }}_2;t_1,t_2\right)=\langle U\left({\mathit{\rho }}_1,t_1\right)U^{*}\left({\mathit{\rho }}_2,t_2\right)\rangle ,$$

(1)

where $U\left(\mathit{\rho },t\right)$ is an instance of a random optical field evaluated at position $\mathit{\rho }=\widehat{\mathit{x}}x+\widehat{\mathit{y}}y$ and time t, * is the complex conjugate, and $\langle \rangle$ is the average over the ensemble of U realizations. The first and second moments of many optical sources vary slowly with time such that they can be considered static. This characteristic is referred to as wide-sense stationary and implies that the MCF defined in Equation (1) depends temporally only on the time difference $\tau =t_1-t_2$.

The temporal Fourier transform of the MCF is the cross-spectral density (CSD) function:

$$W\left({\mathit{\rho }}_1,{\mathit{\rho }}_2;\omega \right)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\Gamma \left({\mathit{\rho }}_1,{\mathit{\rho }}_2;\tau \right)exp\left(\mathrm{j}\omega \tau \right)\mathrm{d}\tau ,$$

(2)

where $\omega$ is the radian frequency. Letting ${\mathit{\rho }}_1={\mathit{\rho }}_2=\mathit{\rho }$ yields the spectral density (SD), namely,

$$S\left(\mathit{\rho };\omega \right)=W\left(\mathit{\rho },\mathit{\rho };\omega \right),$$

(3)

which is the optical intensity at frequency $\omega$. Normalizing the CSD produces the spectral degree of coherence (SDoC), i.e.,

$$\mathsf{\mu }\left({\mathit{\rho }}_1,{\mathit{\rho }}_2;\omega \right)=\frac{W\left({\mathit{\rho }}_1,{\mathit{\rho }}_2;\omega \right)}{\sqrt{S\left({\mathit{\rho }}_1;\omega \right)S\left({\mathit{\rho }}_2;\omega \right)}},$$

(4)

which is a measure of the spatial coherence of the source at frequency $\omega$.

Working with the CSD, as opposed to the MCF, has two significant advantages. The first is the separation of spatial and temporal coherence. In the MCF, time and space are coupled in such a way that determining whether the phenomenon is caused by temporal or spatial coherence can be difficult. On the other hand, the CSD provides a measure of the spatial coherence of an optical source at a certain frequency $\omega$. In this way, the CSD is much more physical than the MCF, and this physical utility has led to the prediction and subsequent experimental validation of many counter-intuitive optical phenomena, most notably, coherence-induced spectral changes [1,2] and coherence-induced polarization changes [3,4,5].

The second advantage is analytical efficiency. The expression for the propagation of mutual coherence is a complicated four-dimensional (4D) integral over space, where space and time are coupled [6]. The typical approach for dealing with this is to assume that the source is narrowband and quasi-monochromatic. These assumptions simplify the mutual coherence propagation expression significantly; however, they restrict the applicability of the results. On the other hand, the expression for the propagation of CSD—in its most general, paraxial form—is a 4D Fresnel transform [1,5]. In many practical scenarios, the 4D Fresnel transform simplifies to a 4D Fourier transform. The results obtained using these propagation integrals are accurate regardless of the source’s bandwidth.

Wolf’s CSD formalism has led to the understanding and prediction of many optical phenomena. Among these, Wolf and others showed that a spatially partially coherent beam (PCB) can be highly directional similar to a laser [1,7,8,9]. Other research soon followed, showing that the reduced spatial coherence of PCBs makes them less susceptible to scintillation or speckle [5,10,11,12,13]. These two discoveries have motivated much of the recent PCS research because high directionality and speckle/scintillation resistance are ideal source characteristics for directed energy, free-space optical communications (FSOC), medical, and numerous other applications.

Because of their many potential uses, the literature is replete with techniques to synthesize sources with controllable spatial coherence. A survey of the literature reveals two general approaches for generating PCBs. The first starts with a spatially coherent source, commonly a laser, and “spoils” the coherence by randomizing the wavefront. This is accomplished using spatial light modulators (SLMs) or rotating ground glass diffusers (GGDs) [5,14,15]. The second approach starts with a spatially incoherent source, an LED for instance, and exploits the Van Cittert–Zernike theorem (VCZT) or uses spatial-frequency filters (typically referred to as just spatial filters) to synthesize the desired source [16,17,18,19,20]. Hybrid approaches have also been developed [11,12,21].

The “coherent” and “incoherent” PCS synthesis methods discussed above have pros and cons that make one more suitable for a specific application than the other. FSOC is a good example and will be the focus of the remainder of this paper because of the author’s research interests. In FSOC, the source’s wavefront cycling rate, i.e., the rate at which the source produces statistically independent wavefronts, must be at least 10–100 times larger than the communication modulation frequency. This is so the detector integrates many independent intensity realizations per digital bit to reduce scintillation caused by atmospheric turbulence between the transmitter and receiver. Considering that FSOC data rates are gigabits to potentially terabits per second and SLM frame rates and GGD rotation rates are optimistically 10 kHz, spoiling the coherence of a laser source is not a realistic option.

On the other hand, the wavefront cycling rate of a spatially incoherent source (e.g., an LED, thermal source, et cetera) is determined by the physics of how the source generates light (for the aforementioned sources, the physical mechanism is spontaneous emission) and is ultimately related to the source’s bandwidth [6]. The bandwidths of such sources applicable to FSOC vary from a couple THz to 10 s of THz. Even for sources at the low end of this range, the wavefront cycles so quickly that 100 gigabit per second data rates are possible. Because of this, the incoherent option is the only viable approach for FSOC applications.

Of the two incoherent techniques, i.e., VCZT and spatial filtering, the latter has many advantages including simplicity and compactness. Surprisingly, despite these advantages and its common use in optical Fourier processing [22], the spatial filtering of a partially coherent source (PCS) has been discussed in relatively few papers [19,20,23,24]. These references, after presenting the requisite Fourier and statistical optics theory, generally focus on specific PCSs, e.g., Gaussian Schell model [1,5] or incoherent sources. Simple, physical, “engineering” expressions for the PCS at the filter output, or directions for designing/choosing a filter to produce a desired beam are not presented.

The goal of this paper is to derive these expressions and demonstrate their use. In the next section, starting with a stochastic source field, the field at the output of a general spatial filter is derived. Taking the autocorrelation of the output field and assuming that the source field can be modeled as a quasi-homogeneous source [1,6] yields closed-form, approximate expressions for the output SD and SDoC. Inverting the SDoC expression produces a function for the filter in terms of the desired output SDoC. Lastly, Section 3 validates these theoretical expressions and demonstrates how to use them to control coherence and beam shape via Monte Carlo wave-optics simulations.

2. Materials and Methods

2.1. Fourier and Statistical Optics Theory

The spatial filter geometry is shown in Figure 1. This simple geometry depicts a random field passing through a two-lens system. The source plane is one focal length $f_1$ in front of ${\mathrm{L}}_1$. The field passes through ${\mathrm{L}}_1$ and then a complex amplitude filter T placed at the rear focus of ${\mathrm{L}}_1$. The location of T also corresponds to the front focal plane of ${\mathrm{L}}_2$. The field then transits ${\mathrm{L}}_2$ before finally being observed in ${\mathrm{L}}_2$’s rear focal plane.

Assuming that ${\mathrm{L}}_1$ is large enough in diameter to collect all the light emitted from the source, the field in ${\mathrm{L}}_1$’s focal plane, right before T, is

$$\begin{array}{cc}\hfill U\left(\mathit{\rho },f_1^{-}\right)& =\frac{exp\left(\mathrm{j}k{f}_1\right)}{\mathrm{j}\lambda f_1}{\int \int }_{-\infty }^{\infty }{U}^{\mathrm{src}}\left({\mathit{\rho }}^{\prime }\right)exp\left(-\mathrm{j}\frac{k}{f_1}\mathit{\rho }·{\mathit{\rho }}^{\prime }\right){\mathrm{d}}^2{\rho }^{\prime }\hfill \\ & =\frac{exp\left(\mathrm{j}k{f}_1\right)}{\mathrm{j}\lambda f_1}{\tilde{U}}^{\mathrm{src}}\left(\frac{\mathit{\rho }}{\lambda f_1}\right)\hfill \end{array},$$

(5)

where $\mathit{\rho }$ is the observation vector, ${\mathit{\rho }}^{\prime }$ is the source vector, $k=2\pi /\lambda$, and $\lambda$ is the wavelength [22]. The random source field $U^{\mathrm{src}}$ is assumed to be a sample function drawn from a wide-sense stationary random process (the dependence of $U^{\mathrm{src}}$ on radian frequency $\omega$ is assumed and suppressed), and ${\tilde{U}}^{\mathrm{src}}$ is the spatial Fourier transform of $U^{\mathrm{src}}$:

$${\tilde{U}}^{\mathrm{src}}\left(\mathit{f}\right)={\int \int }_{-\infty }^{\infty }{U}^{\mathrm{src}}\left(\mathit{\rho }\right)exp\left(-\mathrm{j}2\pi \mathit{f}·\mathit{\rho }\right){\mathrm{d}}^2\rho ,$$

(6)

where $\mathit{f}=\widehat{\mathit{x}}{f}_x+\widehat{\mathit{y}}{f}_y$ is the spatial frequency vector.

The field immediately to the right of the complex amplitude filter is

$$U\left(\mathit{\rho },f_1^{+}\right)=T\left(\mathit{\rho }\right)U\left(\mathit{\rho },f_1^{-}\right)=\frac{exp\left(\mathrm{j}k{f}_1\right)}{\mathrm{j}\lambda f_1}T\left(\mathit{\rho }\right){\tilde{U}}^{\mathrm{src}}\left(\frac{\mathit{\rho }}{\lambda f_1}\right).$$

(7)

Lastly, assuming that ${\mathrm{L}}_2$ is large enough to capture all the light that passes through T, the field at the output of the spatial filter $U^{\mathrm{out}}$ is [22]

$$\begin{array}{cc}\hfill U\left(\mathit{\rho },f_2\right)& =U^{\mathrm{out}}\left(\mathit{\rho }\right)\hfill \\ & =\frac{exp\left(\mathrm{j}k{f}_2\right)}{\mathrm{j}\lambda f_2}{\int \int }_{-\infty }^{\infty }U\left({\mathit{\rho }}^{\prime },f_1^{+}\right)exp\left(-\mathrm{j}\frac{k}{f_2}\mathit{\rho }·{\mathit{\rho }}^{\prime }\right){\mathrm{d}}^2{\rho }^{\prime }\hfill \\ & =-\frac{exp\left[\mathrm{j}k\left(f_1+f_2\right)\right]}{{\lambda }^2{f}_1{f}_2}{\int \int }_{-\infty }^{\infty }T\left({\mathit{\rho }}^{\prime }\right){\tilde{U}}^{\mathrm{src}}\left(\frac{{\mathit{\rho }}^{\prime }}{\lambda f_1}\right)exp\left(-\mathrm{j}\frac{k}{f_2}\mathit{\rho }·{\mathit{\rho }}^{\prime }\right){\mathrm{d}}^2{\rho }^{\prime }\hfill \end{array}.$$

(8)

The output field given in Equation (8) is for a single realization of the random source field $U^{\mathrm{src}}$. Taking the autocorrelation of $U^{\mathrm{out}}$ produces

$$\begin{array}{c}\langle U^{\mathrm{out}}\left({\mathit{\rho }}_1\right)U^{\mathrm{out}*}\left({\mathit{\rho }}_2\right)\rangle =W^{\mathrm{out}}\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)\hfill \\ =\frac{1}{{\left({\lambda }^2{f}_1{f}_2\right)}^2}{\int \int \int \int }_{-\infty }^{\infty }T\left({\mathit{\rho }}_1^{\prime }\right)T^{*}\left({\mathit{\rho }}_2^{\prime }\right){\tilde{W}}^{\mathrm{src}}\left(\frac{{\mathit{\rho }}_1^{\prime }}{\lambda f_1},\frac{{\mathit{\rho }}_2^{\prime }}{\lambda f_1}\right)exp\left[-\mathrm{j}\frac{k}{f_2}\left({\mathit{\rho }}_1·{\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2·{\mathit{\rho }}_2^{\prime }\right)\right]{\mathrm{d}}^2{\rho }_1^{\prime }{\mathrm{d}}^2{\rho }_2^{\prime }\hfill \end{array},$$

(9)

where $W^{\mathrm{out}}$ is the output CSD function [1,5] and ${\tilde{W}}^{\mathrm{src}}$ is the Fourier transform of the source CSD function. It is quite easy to show that an equivalent and more physical expression for $W^{\mathrm{out}}$ is

$$\begin{array}{c}{W}^{\mathrm{out}}\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)={\left[\frac{M}{{\left(\lambda f_2\right)}^2}\right]}^2{\int \int \int \int }_{-\infty }^{\infty }{W}^{\mathrm{src}}\left({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_2^{\prime }\right)\tilde{T}\left(\frac{{\mathit{\rho }}_1-M{\mathit{\rho }}_1^{\prime }}{\lambda f_2}\right){\tilde{T}}^{*}\left(\frac{{\mathit{\rho }}_2-M{\mathit{\rho }}_2^{\prime }}{\lambda f_2}\right){\mathrm{d}}^2{\rho }_1^{\prime }{\mathrm{d}}^2{\rho }_2^{\prime }\hfill \end{array},$$

(10)

where $M=-f_2/f_1$ is the magnification. This expression states that $W^{\mathrm{out}}$ is the scaled convolution of $W^{\mathrm{src}}$ with the Fourier transform of the CSD function of T. In arriving at Equations (9) and (10), no assumptions have been made about the spatial statistics of $U^{\mathrm{src}}$; thus, they are accurate for a source of any state of coherence.

To make further progress, an analytical form for $W^{\mathrm{src}}$ is needed. Because of its use in FSOC research, the physical source modeled here is the output of a multimode fiber (MMF) excited by a “broadband” (1–2 THz) laser [25,26,27]. This source can be approximated by a quasi-homogeneous CSD function [1,5,6], which takes the form

$$W^{\mathrm{src}}\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)=S^{\mathrm{src}}\left(\frac{{\mathit{\rho }}_1+{\mathit{\rho }}_2}{2}\right){\mathsf{\mu }}^{\mathrm{src}}\left({\mathit{\rho }}_1-{\mathit{\rho }}_2\right),$$

(11)

where $S^{\mathrm{src}}$ is the source’s SD, ${\mathsf{\mu }}^{\mathrm{src}}$ is the source’s SDoC, and $S^{\mathrm{src}}$ is a “slow function” (varies slowly) compared to ${\mathsf{\mu }}^{\mathrm{src}}$.

Taking the Fourier transform of Equation (11) and substituting the resulting expression into Equation (9) produces

$$\begin{array}{c}{W}^{\mathrm{out}}\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)=\frac{1}{{\left({\lambda }^2{f}_1{f}_2\right)}^2}{\int \int \int \int }_{-\infty }^{\infty }T\left({\mathit{\rho }}_1^{\prime }\right)T^{*}\left({\mathit{\rho }}_2^{\prime }\right){\tilde{\mathsf{\mu }}}^{\mathrm{src}}\left(\frac{{\mathit{\rho }}_1^{\prime }+{\mathit{\rho }}_2^{\prime }}{2\lambda f_1}\right)\hfill \\ \times {\tilde{S}}^{\mathrm{src}}\left(\frac{{\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime }}{\lambda f_1}\right)exp\left[-\mathrm{j}\frac{k}{f_2}\left({\mathit{\rho }}_1·{\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2·{\mathit{\rho }}_2^{\prime }\right)\right]{\mathrm{d}}^2{\rho }_1^{\prime }{\mathrm{d}}^2{\rho }_2^{\prime }\hfill \end{array}.$$

(12)

Making the variable substitutions $\mathit{s}={\mathit{\rho }}_1^{\prime }-{\mathit{\rho }}_2^{\prime }$ and $\mathit{t}=\left({\mathit{\rho }}_1^{\prime }+{\mathit{\rho }}_2^{\prime }\right)/2$ and simplifying, yields

$$\begin{array}{c}{W}^{\mathrm{out}}\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)=\frac{1}{{\left({\lambda }^2{f}_1{f}_2\right)}^2}{\int \int }_{-\infty }^{\infty }{\tilde{\mathsf{\mu }}}^{\mathrm{src}}\left(\frac{\mathit{t}}{\lambda f_1}\right)exp\left[-\mathrm{j}\frac{k}{f_2}\left({\mathit{\rho }}_1-{\mathit{\rho }}_2\right)·\mathit{t}\right]\hfill \\ \times {\int \int }_{-\infty }^{\infty }T\left(\mathit{t}+\frac{1}{2}\mathit{s}\right)T^{*}\left(\mathit{t}-\frac{1}{2}\mathit{s}\right){\tilde{S}}^{\mathrm{src}}\left(\frac{\mathit{s}}{\lambda f_1}\right)exp\left[-\mathrm{j}\frac{k}{f_2}\left(\frac{{\mathit{\rho }}_1+{\mathit{\rho }}_2}{2}\right)·\mathit{s}\right]{\mathrm{d}}^{2}s{\mathrm{d}}^{2}t\hfill \end{array}.$$

(13)

Now, assuming ${\tilde{S}}^{\mathrm{src}}$ is a very narrow function of $\mathit{s}$ (much narrower than T), Equation (13) can be approximated as

$$\begin{array}{c}{W}^{\mathrm{out}}\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)\approx \frac{1}{{\left({\lambda }^2{f}_1{f}_2\right)}^2}{\int \int }_{-\infty }^{\infty }{\left|T\left(\mathit{t}\right)\right|}^2{\tilde{\mathsf{\mu }}}^{\mathrm{src}}\left(\frac{\mathit{t}}{\lambda f_1}\right)exp\left[-\mathrm{j}\frac{k}{f_2}\left({\mathit{\rho }}_1-{\mathit{\rho }}_2\right)·\mathit{t}\right]{\mathrm{d}}^{2}t\hfill \\ \times {\int \int }_{-\infty }^{\infty }{\tilde{S}}^{\mathrm{src}}\left(\frac{\mathit{s}}{\lambda f_1}\right)exp\left[-\mathrm{j}\frac{k}{f_2}\left(\frac{{\mathit{\rho }}_1+{\mathit{\rho }}_2}{2}\right)·\mathit{s}\right]{\mathrm{d}}^{2}s\hfill \end{array}.$$

(14)

Note that this assumption is equivalent to assuming that the integrand of Equation (12), without the Fourier kernel, is another quasi-homogeneous source. Evaluating the above Fourier integrals produces

$$\begin{array}{c}{W}^{\mathrm{out}}\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)=\frac{1}{M^2}{S}^{\mathrm{src}}\left(\frac{{\mathit{\rho }}_1+{\mathit{\rho }}_2}{2M}\right){\mu }^{\mathrm{filt}}\left({\mathit{\rho }}_1-{\mathit{\rho }}_2\right)\hfill \\ {\mathsf{\mu }}^{\mathrm{filt}}\left({\mathit{\rho }}_1-{\mathit{\rho }}_2\right)={\int \int }_{-\infty }^{\infty }{\left|T\left(\lambda f_1\mathit{u}\right)\right|}^2{\tilde{\mathsf{\mu }}}^{\mathrm{src}}\left(\mathit{u}\right)exp\left[\mathrm{j}2\pi \left(\frac{{\mathit{\rho }}_1-{\mathit{\rho }}_2}{M}\right)·\mathit{u}\right]{\mathrm{d}}^{2}u\hfill \end{array}.$$

(15)

The output SD and SDoC are

$$\begin{array}{c}{S}^{\mathrm{out}}\left(\mathit{\rho }\right)=W^{\mathrm{out}}\left(\mathit{\rho },\mathit{\rho }\right)=\frac{1}{M^2}{S}^{\mathrm{src}}\left(\frac{\mathit{\rho }}{M}\right){\mu }^{\mathrm{filt}}\left(0\right)\hfill \\ {\mathsf{\mu }}^{\mathrm{out}}\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)=\frac{W^{\mathrm{out}}\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)}{\sqrt{S^{\mathrm{out}}\left({\mathit{\rho }}_1\right)S^{\mathrm{out}}\left({\mathit{\rho }}_2\right)}}\approx \frac{{\mathsf{\mu }}^{\mathrm{filt}}\left({\mathit{\rho }}_1-{\mathit{\rho }}_2\right)}{{\mathsf{\mu }}^{\mathrm{filt}}\left(0\right)}\hfill \end{array}.$$

(16)

Note that the output SD $S^{\mathrm{out}}$ (physically, the source shape) is a magnified version of $S^{\mathrm{src}}$, and the output SDoC ${\mathsf{\mu }}^{\mathrm{out}}$ is a magnified and filtered version of ${\mathsf{\mu }}^{\mathrm{src}}$.

Although expressions for $S^{\mathrm{out}}$ and ${\mathsf{\mu }}^{\mathrm{out}}$ have been derived, an equation for T in terms of ${\mathsf{\mu }}^{\mathrm{out}}$ (i.e., what T yields a given ${\mathsf{\mu }}^{\mathrm{out}}$) is most useful. Starting with the above expression for ${\mathsf{\mu }}^{\mathrm{out}}$ given in Equation (16) and assuming that ${\tilde{\mathsf{\mu }}}^{\mathrm{src}}$ is a slow function compared to T (this is actually a consequence of $W^{\mathrm{src}}$ assuming a quasi-homogeneous form) produces

$${\mathsf{\mu }}^{\mathrm{out}}\left({\mathit{\rho }}_1-{\mathit{\rho }}_2\right)\approx \frac{{\int \int }_{-\infty }^{\infty }{\left|T\left(\lambda f_1\mathit{u}\right)\right|}^{2}exp\left[\mathrm{j}2\pi \left(\frac{{\mathit{\rho }}_1-{\mathit{\rho }}_2}{M}\right)·\mathit{u}\right]{\mathrm{d}}^{2}u}{{\int \int }_{-\infty }^{\infty }{\left|T\left(\lambda f_1\mathit{u}\right)\right|}^2{\mathrm{d}}^{2}u}.$$

(17)

Letting $\mathbf{\Delta }\mathit{\rho }={\mathit{\rho }}_1-{\mathit{\rho }}_2$ for convenience and substituting $\mathit{t}=\lambda f_1\mathit{u}$ into both the numerator and denominator integrals simplifies Equation (17) to

$${\mathsf{\mu }}^{\mathrm{out}}\left(\mathbf{\Delta }\mathit{\rho }\right)=\frac{{\int \int }_{-\infty }^{\infty }{\left|T\left(\mathit{t}\right)\right|}^{2}exp\left[-\mathrm{j}2\pi \left(\frac{\mathbf{\Delta }\mathit{\rho }}{\lambda f_2}\right)·\mathit{t}\right]{\mathrm{d}}^{2}t}{{\int \int }_{-\infty }^{\infty }{\left|T\left(\mathit{t}\right)\right|}^2{\mathrm{d}}^{2}t}.$$

(18)

Applying the Fourier transform defined in Equation (6) to both sides of Equation (18) produces the desired result:

$$\frac{{\left|T\left(-\lambda f_2\mathit{f}\right)\right|}^2}{{\int \int }_{-\infty }^{\infty }{\left|T\left(\lambda f_2\mathit{f}\right)\right|}^2{\mathrm{d}}^{2}f}={\int \int }_{-\infty }^{\infty }{\mathsf{\mu }}^{\mathrm{out}}\left(\mathbf{\Delta }\mathit{\rho }\right)exp\left(-\mathrm{j}2\pi \mathit{f}·\mathbf{\Delta }\mathit{\rho }\right){\mathrm{d}}^2\Delta \rho ={\mathsf{\Phi }}^{\mathrm{out}}\left(\mathit{f}\right).$$

(19)

This expression states that the spatial power spectrum of the output field ${\mathsf{\Phi }}^{\mathrm{out}}$ (the Fourier transform of ${\mathsf{\mu }}^{\mathrm{out}}$) is proportional to the magnitude squared of T rotated 180${}^{\circ }$. The proportionality constant is there to normalize the area under the power spectrum; it is inconsequential in practice. The above interpretation of Equation (19) leads to the following steps for finding T:

• Choose the desired output SDoC ${\mathsf{\mu }}^{\mathrm{out}}$.
• Compute the spatial Fourier transform of ${\mathsf{\mu }}^{\mathrm{out}}$ to find ${\mathsf{\Phi }}^{\mathrm{out}}$.
• Take the square root of ${\mathsf{\Phi }}^{\mathrm{out}}$ to find the magnitude of T.
• Normalize T such that $\underset{\mathit{\rho }\in {\mathbb{R}}^2}{max}\left|T\left(\mathit{\rho }\right)\right|=1$.
• Rotate T 180${}^{\circ }$.

It should be stated that Equations (16) and (19) are accurate if the following three conditions hold:

• The diameters of lenses ${\mathrm{L}}_1$ and ${\mathrm{L}}_2$ are large enough to collect all the light emitted from the source and the light that passes through the spatial filter T, respectively.
• The random source field $U^{\mathrm{src}}$ is well approximated by a CSD function with a quasi-homogeneous form.
• ${\tilde{S}}^{\mathrm{src}}$ is a fast function (narrow function) compared to T.

The first item is easily met in practice by using ray tracing and other standard optical design techniques. The second criterion is true for “spatially incoherent” sources such as LEDs or, applicable here, the light emitted from a broadband-laser-excited MMF. The third condition is the hardest to satisfy in practice. The effects of violating this condition are shown in Section 3.1.

Because of the assumptions that underpin the analysis (enumerated above), Equations (16) and (19) are approximations of the output PCS moments. Recall that the goal here was to find simple, physical, “engineering” expressions that could be used to design a spatial filter that would produce a desired PCS. Equations (16) and (19) achieve that objective.

2.2. Simulation Details

The purpose of the Monte Carlo wave-optics simulations is to validate the above analysis and to demonstrate how to apply it, specifically Equations (16) and (19), to generate a desired PCS. Before proceeding to the results, a brief discussion of the set-up is warranted.

To simulate the output of a MMF fed by a broadband laser [25,26,27], $W^{\mathrm{src}}$ took the form

$$W^{\mathrm{src}}\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)=circ\left(\frac{{\rho }_1}{D_c/2}\right)circ\left(\frac{{\rho }_2}{D_c/2}\right)jinc\left(\frac{k}{2}\mathrm{NA}\left|{\mathit{\rho }}_1-{\mathit{\rho }}_2\right|\right),$$

(20)

where $circ\left(x\right)$ is the circle function defined by Goodman [22], $jinc\left(x\right)=2J_1\left(x\right)/x$, and $J_1$ is a first-order Bessel function of the first kind. The fiber’s core diameter and numerical aperture were $D_c=105\mathsf{\mu }\mathrm{m}$ and $\mathrm{NA}=0.22$, respectively. The simulated wavelength was $\lambda =635\mathrm{nm}$. Instances of $U^{\mathrm{src}}$ were synthesized from $W^{\mathrm{src}}$ using the method described in [28].

The focal lengths of the lenses comprising the spatial filter were $f_1=100\mathrm{mm}$ and $f_2=500\mathrm{mm}$ (see Figure 1). The source, focal, and output planes were discretized using 512 points per side with grid spacings equal to 0.8203 $\mathsf{\mu }$m, 0.1512 mm, and 4.102 $\mathsf{\mu }$m, respectively. All propagations were performed using fast Fourier transforms [29].

The desired, output SDoC ${\mathsf{\mu }}^{\mathrm{out}}$ (after filtering) was

$${\mathsf{\mu }}^{\mathrm{out}}\left({\mathit{\rho }}_1-{\mathit{\rho }}_2\right)={sinc}^2\left(\pi \frac{x_1-x_2}{2D_c}\right){sinc}^2\left(\pi \frac{y_1-y_2}{4D_c}\right),$$

(21)

where $sinc\left(x\right)=sin\left(x\right)/x$. This SDoC was chosen to show that a source with a rotationally invariant CSD function (see Equation (20)) can be transformed (via filtering) into one with different, specified x and y correlation lengths. The required T to synthesize a source with this SDoC was determined using the five steps enumerated above. Figure 2a shows this T.

The simulated output SD $S^{\mathrm{out},\mathrm{sim}}$ and ${\mathsf{\mu }}^{\mathrm{out},\mathrm{sim}}$ were computed from 2000 propagated instances of $U^{\mathrm{src}}$. These results were then compared to the theoretical predictions given in Equation (16). All simulations were performed using MATLAB^® version R2017a (The MathWorks, Inc., Natick, MA, USA); the scripts (.m files) are included as Supplementary Materials.

3. Results and Discussion

3.1. Spatial Coherence Control

Figure 3 shows the simulation results. Figure 3a,b shows the theoretical and simulated SDs ($S^{\mathrm{out},\mathrm{thy}}$ obtained from Equation (16) and $S^{\mathrm{out},\mathrm{sim}}$); Figure 3c,d shows the theoretical and simulated SDoCs (${\mathsf{\mu }}^{\mathrm{out},\mathrm{thy}}$ given in Equation (21) and ${\mathsf{\mu }}^{\mathrm{out},\mathrm{sim}}$); Figure 3e,f shows the $y=0$ and $x=0$ slices through $S^{\mathrm{out},\mathrm{thy}}$ and $S^{\mathrm{out},\mathrm{sim}}$; and Figure 3g,h shows the $y=0$ and $x=0$ slices through ${\mathsf{\mu }}^{\mathrm{out},\mathrm{thy}}$ and ${\mathsf{\mu }}^{\mathrm{out},\mathrm{sim}}$.

The $S^{\mathrm{out},\mathrm{thy}}$ and $S^{\mathrm{out},\mathrm{sim}}$ results noticeably disagree. This occurs because ${\tilde{S}}^{\mathrm{src}}$ is not a fast function compared to T, viz., Criterion 3 is violated. Recall that this assumption was made in going from Equation (13) to Equation (14). The width of ${\mathsf{\mu }}^{\mathrm{out}}$ in the y direction is larger than in the x direction (see Equation (21)). The Fourier transform relationship between T and ${\mathsf{\mu }}^{\mathrm{out}}$ means T has complementary widths to ${\mathsf{\mu }}^{\mathrm{out}}$, i.e., T is wider in the x direction than in the y direction. This ${\mathsf{\mu }}^{\mathrm{out}}$-T relationship explains why $S^{\mathrm{out},\mathrm{sim}}$ compares better to $S^{\mathrm{out},\mathrm{thy}}$ in Figure 3e than in Figure 3f.

Although $S^{\mathrm{out},\mathrm{thy}}$ and $S^{\mathrm{out},\mathrm{sim}}$ disagree, the objective here is to control spatial coherence as measured by the SDoC. Figure 3c,d,g,h shows these results. Clearly, the effects of violating Criterion 3 are generally limited to the output SD because the agreement between ${\mathsf{\mu }}^{\mathrm{out},\mathrm{thy}}$ and ${\mathsf{\mu }}^{\mathrm{out},\mathrm{sim}}$ is excellent. These results validate the PCS spatial filtering method presented here.

3.2. Beam Shaping

The results presented in Figure 3 demonstrated that the output SDoC could be precisely controlled by the proper choice of T. A related application where the Section 2.1 analysis can also be applied is beam shaping.

The goal in beam shaping is to produce a beam with a desired SD at a desired location—typically, at the focus of a lens or, equivalently, in the far zone of the source. Referring back to Equation (15), propagating $W^{\mathrm{out}}$ to the far field, and setting ${\mathit{\rho }}_1={\mathit{\rho }}_2=\mathit{\rho }$ yields the far-zone SD:

$$S^{\mathrm{fz}}\left(\mathit{\rho },z\right)\approx \frac{M^2}{{\lambda }^2{z}^2}{\tilde{S}}^{\mathrm{src}}\left(0\right){\tilde{\mathsf{\mu }}}^{\mathrm{src}}\left(0\right){\left|T\left(-f_2\frac{\mathit{\rho }}{z}\right)\right|}^2.$$

(22)

Note that the source takes the shape of the absolute square of the filter—magnified and rotated 180${}^{\circ }$—in the far field. This can be exploited to produce beams of any desired shape.

As an example, Figure 4 shows the results of a beam shaping simulation. The source (see Equation (20)) and set-up were identical to that described in Section 2.2, only T was changed to produce a beam with a more complex shape. Figure 2b shows this T.

Figure 4a,b shows $S^{\mathrm{out},\mathrm{thy}}$ (obtained from Equation (16)) and $S^{\mathrm{out},\mathrm{sim}}$, respectively. By design, these SDs should look very similar to those in Figure 3a,b—the differences are due to the different T. Figure 4c shows the absolute value of the pixel-by-pixel difference of Figure 4a,b.

Figure 4d,e shows the theoretical and simulated far-zone SDs, $S^{\mathrm{fz},\mathrm{thy}}$ obtained from Equation (22) and $S^{\mathrm{fz},\mathrm{sim}}$, respectively. The theoretical far-zone SD $S^{\mathrm{fz},\mathrm{thy}}$ is the desired beam shape. $S^{\mathrm{fz},\mathrm{sim}}$ was computed from 2000 far-zone propagated instances of $U^{\mathrm{out}}$ (these propagations were simulated using fast Fourier transforms [29]), which were obtained from 2000 instances of $U^{\mathrm{src}}$, synthesized in the manner described in [28]. Similarly, Figure 4c,f shows the absolute value of the pixel-by-pixel difference of Figure 4d,e.

The agreement between the simulated and theoretical SDs is excellent. These results further validate the analysis presented in Section 2.1.

If one were to physically generate this source, it would start as partially coherent light emitted from a MMF. It would exit the spatial filter as a circular flat-topped beam, as shown in Figure 4a,b. As the beam propagated away from the spatial filter output plane, it would morph into the Celtic cross image shown in Figure 4d,e.

4. Conclusions

The control of an optical source’s SDoC via spatial filtering was demonstrated. Starting with a quasi-homogeneous random field, the optical field at the output of a general spatial filter was derived. Taking moments of the stochastic output field yielded closed-form expressions for the output SD and SDoC. The key relation which resulted from this analysis was a closed-form expression for the filter function T in terms of the desired SDoC.

The analytical SD and SDoC equations derived in Section 2.1 were verified by Monte Carlo analysis. First, the T required to produce a specified SDoC was found using the steps developed and discussed in Section 2.1. Then, simulating a source used in FSOC research (namely, the stochastic field emitted from a MMF excited by a broadband laser source) and a simple spatial filter containing T, the output SD and SDoC were computed from many independent realizations of the source field. Lastly, the simulated output SD and SDoC were compared to the theoretical predictions. The agreement between the two was generally good with some discrepancies between the simulated and theoretical SDs. The reason for these differences was discussed.

In addition, a beam shaping simulation was also performed to demonstrate another use of coherence control via spatial filtering. Using the analytical relations derived in Section 2.1, an approximate expression for the desired far-zone SD, in terms of the filter function T, was derived. This equation was used to find the T required to synthesize a beam with a complex shape. The simulated results were in excellent agreement with the derived, analytical SDs.

The spatial filtering technique presented in this paper has many advantages over other PCS synthesis methods. First and foremost among them is speed. As discussed in the Introduction, the wavefront cycling rate for this method can be 10s of THz, which is 10${}^9$ times faster than SLM- or GGD-based approaches. When compared against other incoherent synthesis techniques (e.g., VCZT), spatial filtering is simpler to implement. As demonstrated in Section 3, changing the output field SDoC requires only a filter change. These pros make PCS spatial filtering useful in any application where coherence control is advantageous. These include, but are not limited to, FSOC, directed energy, medicine, manufacturing, and remote sensing.