Synthesizing Time-evolving Partially-coherent Schell-model Sources

Time-evolving simulation of sources with partial spatial and temporal coherence is sometimes instructive or necessary to explain optical coherence effects. Yet, existing time-evolving synthesis techniques often require prohibitive amounts of computer memory. This paper discusses three methods for the synthesis of continuous or pulsed time-evolving sources with nearly arbitrary spatial and temporal coherence. One method greatly reduces computer memory requirements, making this type of synthesis more practical. The utility of all three methods is demonstrated via a modified form of Young's experiment. Numerical simulation and laboratory results for time-averaged irradiance are presented and compared with theory to validate the synthesis techniques.

Considering their many possible applications, techniques to physically realize partially-coherent sources with desired coherence properties are quite numerous [27,8,14,15,[28][29][30][31]9,[32][33][34][35][36]. One approach uses the coherent mode representation to generate ensembles of random pulses, where each pulse is independent of all other pulses [35,36]. More commonly, partially-coherent sources are produced by synthesizing a set of independent, spatially-correlated complex transmittance or traditional phase screens. The desired moments (often the mean irradiance or complex degree of coherence) are computed over the set of random field instances, thus approximating the ensemble statistics. Assuming ergodicity, this approach produces the long-time field moments, which are sufficient in many scenarios. However, in some cases, modeling the time evolution of the partially-coherent source is important. Examples include high-speed detection of shorttime interference, propagation of optical vortices through turbulence, and propagation of excimer beams through optical systems [37][38][39]. In addition, visualizing the time evolution of the source can yield physical insight not provided by the ensemble averaging techniques cited above.
Compared to the ensemble averaging approaches, techniques that model the temporal evolution of partially-coherent sources are relatively few [39,37,40,38]. Of those cited, the works of Rydberg and Bengtsson [40] and Davis [38] are the most relevant to this paper. Both the method of Rydberg and Bengtsson and "the separable case" in Davis produce an instance of a field with a customizable mutual coherence function (MCF) by filtering a three-dimensional (3D) array of delta-correlated circular complex Gaussian random numbers. For computational efficiency, both groups of researchers perform the filtering operation in the frequency domain using the fast Fourier transform (FFT) algorithm. The approach discussed in Refs. [40,38] is very intuitive and permits a broad range of spatial and temporal coherence functions, but it requires a large amount of computer resources (in particular, memory) [40]. Indeed, although Rydberg and Bengtsson do present some results using the approach, the authors generally dismiss it as being too costly [40].
In this paper, three techniques to produce time-evolving partiallycoherent sources are presented. The sources need not be crossspectrally pure or stationary. The first method, the complex screen approach, is essentially equivalent to the technique discussed in Refs. [40,38]. It is presented here for completeness, and more importantly, because the details are necessary in the development of the other two techniques, which address two of the major shortcomings of the complex screen approach. The first of these latter two methods is a phase-only (i.e., a traditional phase screen) technique. The phase screen method, easily implemented using a phase-only spatial light modulator (SLM), is well suited for applications in which control of field amplitude (required in the complex screen approach) is either difficult or undesirable. While not new [8], it has not been previously demonstrated for time-evolving source synthesis, and it is needed for the experimental portion of this work. The final approach is a novel hybrid FFT-convolution variant which is applicable to both the complex screen and phase screen techniques. The hybrid FFT-convolution method significantly reduces the computer memory required to perform a time-evolving partially-coherent source simulation or experiment while maintaining good numerical efficiency. For certain problems, this method avoids the need for resource-intensive distributedmemory supercomputers, thus offering significant practical benefits.
Section 2 presents the theory necessary to implement all three techniques for synthesizing time-evolving partially-coherent fields. Section 3 then presents numerical and laboratory results from a modified Young's interference experiment. Finally, the time-averaged irradiance patterns are compared to theory for validation.

Coherence functions
To begin the derivation, consider a scalar source field , v is the mean frequency, and U is assumed to vary slowly relative to the complex carrier πvt exp(−j2 ). Propagation is in the z direction. The slowly-varying field U is a complex envelope function which defines the spatial and temporal coherence of the source. Thus, the proceeding analysis considers only U. The MCF of U takes the following Schell-model form where γ is the complex degree of coherence (DoC), I is the average irradiance, τ t t = − 1 2 , and ρ ρ ρ Δ = − 1 2 . Before discussing the hybrid approach, we present two methods for generating time-evolving partially-coherent sources. The first method, called the complex screen (CS) approach, assumes a field of the form where T is an instance of a random field which varies in both amplitude and phase. The second method, the phase screen (PS) approach, only requires control of the field phase. The associated field is where ϕ is the phase of the random field instance. Taking the cross-correlations of Eqs.
The angle brackets 〈〉 represent the averaging operator. Because we are working with the complex field envelope U, the above expressions are centered about zero temporal frequency rather than the mean optical frequency. Later, we will observe that the present treatment is justified and greatly eases the numerical analysis. Note that for the complex screen approach, the cross-correlation of T directly relates to γ CS . Thus, to synthesize partially-coherent fields using the complex screen approach, one only needs to generate random numbers whose correlation matches γ CS . This is shown in the next section. On the other hand, for the phase screen approach, the crosscorrelation function of ϕ is not directly related to γ PS . Thus, the derivation of the phase screen approach is more involved than the complex screen method.
As an example of the extra steps required for the phase screen method, first assume the field is cross-spectrally pure such that γ PS is separable in space and time. The phase cross-correlation function γ ϕ PS must also be separable. Let the spatial correlation of both the phase and the field be Gaussian; let the temporal correlation be Lorentzian. The DoC of the field can then be written as where δ is the spatial coherence radius and v Δ is the half-power bandwidth. For the phase ϕ, let the cross-correlation function be where l ϕ and τ ϕ are, in essence, the phase spatial coherence radius and coherence time, respectively. We desire to relate γ PS to γ ϕ PS . As such, assume that the phase is Gaussian distributed with zero mean. Then, Eq. (6) simplifies to where σ ϕ 2 is the phase variance [8,41].
To progress, let σ ⪢1 ) if γ PS is to have significant value. This motivates expanding γ ϕ PS in a Taylor series about zero. Keeping only the first two terms of the expansion yields Now, γ PS simplifies to the desired product of Gaussian and Lorentzian functions: By comparing Eqs. (6) and (10), γ PS achieves the desired form if and Note that the relations are proportional to the phase variance σ ϕ 2 or its square root. So, l ϕ and τ ϕ increase with σ ϕ and σ ϕ 2 , respectively, increasing the width of γ ϕ PS . Thus, for a given spatial coherence radius δ and half-power bandwidth v Δ , the grid size required to adequately sample γ ϕ PS without "clipping" will become excessive as σ ϕ 2 grows.
However, a large value for σ ϕ 2 is needed to achieve good accuracy with the truncated Taylor series approximation. The authors propose σ π = 4 ϕ 2 2 as a good general compromise. For this σ ϕ 2 , the root-meansquare difference between the true γ PS [Eq. (8)] and the approximate γ PS [Eq. (10)] is only about 0.005. From this example, the phase screen method clearly allows exponential forms for γ PS , e.g., Gaussian and Lorentzian. In fact, the exponential is limited to powers less than two, as higher powers yield power spectral densities (PSDs) with non-physical negative values. Of note, Wang et al. showed that other (non-exponential) coherence function forms are possible using Gaussian phase screens [42]. However, this flexibility comes at the expense of light [43].

Source generation
Let us now consider source synthesis starting with the complex screen method. We will compare the complex degree of coherence to its discrete equivalent. By taking the Fourier series expansion of T, which is now assumed to be periodic, where mni are the Fourier series coefficients, and L x , L y , and L t are the periods in x, y, and t, respectively. The DoC of this expansion is Also consider the forward Fourier transform of γ CS , where f is the spatial frequency vector and v is the temporal frequency. This expression is the spatial and temporal PSD, which must be both real and positive. The inverse relation is Eq. (16) must equal Eq. (14). Equating the two [after discretizing Eq. (16)] reveals the autocorrelation of the Fourier series coefficients in terms of the PSD of the field, where n L = / y y , and δ is the discrete Dirac delta function. Thus, the variance of the Fourier series coefficients is The generating function of T is then where x, y, and t are discretized as x gL N = / x x , y hL N = / y y , and t lL N = / t t , and r is a delta-correlated matrix of circular complex Gaussian random numbers with zero mean and unit variance in each of the real and imaginary parts. Note that Eq. (19) is a threedimensional inverse discrete Fourier transform. So, we can use the numerically efficient FFT algorithm to generate T. Alternately, by using the convolution theorem, we can write Eq. (19) as a convolution operation in the physical domain. Thus, the inverse Fourier transform of the positive square root of the PSD acts as a linear shift-invariant filter [40,38]. The filtering operation preserves the Gaussian distribution of r so that the source is thermal or chaotic light, such as a highly multi-mode laser [44]. The total field formed from T is given by Eq. (2). In this equation, ρ I t ( , ) is the beam shape in space and time, e.g., a pulsed Gaussian beam.
Field synthesis using the phase screen method is very similar. In Section 2.1, we related γ PS to γ ϕ PS . Now, following the same general procedure used in the derivation of the complex screen method, we write the Fourier series expansion of the phase and compute its crosscorrelation. Let the PSD of the phase be Φ ϕ PS . Equating the Fourier series expansion of γ ϕ PS to γ ϕ PS , written as the inverse Fourier transform of Φ ϕ PS , we solve for the variance of the Fourier series coefficients. Using this variance in the Fourier series expansion produces Note that only the real part of Eq. (20) is retained here. This is because ϕ must be real. One could also use the imaginary part of Eq. (20). Because r is circular complex Gaussian, both the real and imaginary parts are independent and identically distributed. Like Eqs. (19), (20) can be efficiently evaluated using the FFT. The time-evolving source produced by this phase field is given in Eq. (3).

Hybrid FFT-convolution
Thus far, we have presented two methods for synthesizing timeevolving partially-coherent sources. They involve spectral filtering operations, allowing use of the efficient FFT algorithm. Because T is a 3D matrix, the memory requirements are considerable for wide sources with many time samples. The third method we present is a hybrid FFT-convolution approach. The 2D spatial filtering is still applied using the FFT, but the 1D time filtering is applied as a convolution operation which typically only requires the screens from tens of previous time steps. Thus, a much smaller 3D matrix is stored. At each simulated time step, the matrix is shifted, the oldest screen is removed, and a newly generated screen replaces it, as illustrated in Fig. 1. This hybrid FFT-convolution approach can be applied to both the complex screen and phase screen methods. Because the convolution operation is O N ( ) 2 for an array with N elements, while the FFT is O N N ( log ), this hybrid approach is less efficient than full spectral filtering. Thus, memory requirements are greatly reduced at the expense of computation time.
The convolution filter's scaling and coefficients were derived previously [45]. For brevity, the derivation is not shown here. Consider an instance of a discrete random field T formed by the time convolution of an instance of a unit-variance, zero-mean random field r with a filter function W, where r is spectrally filtered in space and delta-correlated in time. The time spacing is t Δ . Assume that the temporal PSD t is separable from the spatial PSD-an oft-used assumption which does limit this method to cross-spectrally pure fields. The filter function needed to realize the desired t can be found by manipulating the autocorrelation of T, resulting in where −1 represents the inverse Fourier transform. This discrete filter produces a time autocorrelation consistent with t . Because of the Fourier transform relationship between the convolution filter and the PSD, the width of the filter is on the order of one coherence time. So, the stored field's number of time steps is comparable to the number of samples per coherence time. The exact relation depends on the source's lineshape. Because most simulations run over many coherence times, this method greatly reduces the size (number of time steps) of the field stored to memory. The hybrid-FFT convolution process is summarized as follows: 1. Initialize the inputs to the time filtering process by applying spectral filtering in space to the required number of 2D matrices of deltacorrelated circular complex Gaussian random numbers ρ r t ( , ) i . 2. For the current time step, generate one new spatial screen through spectral filtering.
• Shift the new screen into memory and remove the oldest screen.
3. Compute T using the discrete convolution expression given in Eq. (21). 4. Repeat steps 2-4 once for each time step of the simulation.

Young's experiment
To demonstrate the scope and utility of the time-evolving source generation methods presented here, we use a modified version of Young's classic interference experiment, both for numerical simulations and a laboratory experiment. The layout is shown in Fig. 2. In this section, an analytical result for the propagation of synthesized timeevolving fields through this layout is compared to a rigorously-derived solution to validate our treatment of the source.
Consider the geometry depicted in Fig. 2. The light from the source is truncated by two circular pinholes of diameter d and separation D in an opaque screen. The path to one pinhole is optically longer than the other path by distance Δ. This optical path difference (OPD) allows temporal coherence to affect the interference pattern, even for points along the optical axis. The light transmitted by the two pinholes is focused onto the observation plane by an achromatic lens of focal length f.
By rigorously applying optical coherence theory, the time-averaged observation plane irradiance is [2,44] where x J x x jinc( ) = 2 ( )/ 1 [46], I x y ( , ) is the mean irradiance of the source field, t c = Δ/ d is the time delay of the longer path, and c is the speed of light in vacuum. Here, the field is assumed to be wide-sense stationary, Schell-model, and quasimonochromatic. Also, the spatial coherence diameter is assumed to be much larger than the pinholes.
Unlike this rigorous solution, the derivations of the source generation methods considered only the complex envelope, not the total optical field. Deriving the average irradiance from the complex envelope changes the result slightly. In this case, the field just past the pinhole mask is where circ is the circle function with radius d /2 [46]. Recall that U is the complex envelope, which is centered about zero temporal frequency. Because of this, the factor πvt exp(j2 ) d is artificially added. This factor adds the proper "fringe steering" to the final result while avoiding the need to sample the optical frequency, which, because of the Nyquist sampling criterion and v v ⪢Δ , is numerically unpalatable. The timeaveraged observation plane irradiance 〈I 〉 sim can now be found by propagating U in Eq. (24) to the focal/observation plane, taking the magnitude squared of the resulting field to derive the instantaneous irradiance, and then time averaging. The resulting expression is   Young's interference experiment diagram. A source of arbitrary spatial and temporal coherence is truncated by an opaque mask with two circular pinholes. The light reaching one pinhole experiences an optical path difference Δ. An achromatic lens focuses all wavelengths toward the same point in the observation plane. Note that this interference experiment was also analyzed in Ref. [40].

Optics Communications 387 (2017) 377-384
Note that there is only one difference between I theory and I sim . Namely, γ Re[ ] in I theory includes a Dx fc /( ) term, while the same factor in I sim does not. This term accounts for the possible reduction in coherence due to the OPD between the paths from each pinhole and the observation point. The loss of this Dx fc /( ) term in I sim is caused by neglecting the mean optical frequency.
At this point, it is logical to consider whether neglecting Dx fc /( ) is acceptable. The irradiance pattern I theory is the classic Airy pattern plus fringes due to interference. Because the side lobes of the Airy pattern are much weaker than the main lobe, a vast majority of the energy is contained within the width of the main jinc lobe. Thus, the effective maximum value of Dx fc /( ) is that term evaluated at λf d 1.

Numerical young's experiment
Let us now consider the numerical propagation of source fields through the modified Young's experiment shown in Fig. 2. The source is generated using the complex screen method. The time delay is realized by referencing an earlier time step of the source for the delayed path compared to that used for the non-delayed path. To propagate the field from the source plane to the observation plane, the Fresnel diffraction integral is discretized and executed using the FFT algorithm in a computer. This process is repeated for each time step of the simulation. Thus, the observation plane irradiance pattern is produced as a function of time. Taking the average over many coherence times produces an irradiance pattern which can be compared to theory [Eq. (23)]. Fig. 3 shows the comparisons of theory and the time-averaged simulation results. Here, cross-section plots of the irradiance are shown along x at y=0. The conditions are uniform 1 W/m 2 mean source irradiance, 1 μm wavelength, 256 grid points per side, 3 cm spatial grid extent, and 30,000 time steps over 100 ns. Cross spectral purity is assumed. The line shape of the source is rectangular with v Δ = 30GHz; the associated coherence time is 33.3 ps [44]. The source's spatial coherence function is Gaussian with a δ e = 0.707 mm1/ (coherence) radius. The DoC of the source is where x π x π x sinc( ) = sin( )/( ) [46]. The pinholes are 117 μm in diameter and are separated by 352 μm. The lens focal length is 1 m. From left to right in Fig. 3, the time delays t d are 16.7 ps, 33.3 ps, and 50.0 ps, respectively.
Theory and simulation results agree very well. The overall patterns are Airy disks due to the circular shape of the pinholes. Within the Airy disks, fringes of constructive and destructive interference form. Complete nulls will not occur, even if Δ = 0, as D δ ≈ ; thus, spatial coherence reduces fringe visibility. In all cases shown, t d is comparable to the coherence time, so fringe visibility is further reduced by temporal coherence. Note that fringe visibility vanishes for the 33.3 ps delay (Fig. 3b), which equals one coherence time. This value is the first null of the sinc temporal coherence function. Fringes reappear in Fig. 3c, as the time delay falls between the first and second nulls of the sinc function. Also note that the fringe locations change with the time delay. This steering in the numerical results matches the theory quite well due to the inclusion of πvt exp(j2 ) d in Eq. (24). In Fig. 3a, numerical complex screen results are shown for both the standard implementation and the hybrid FFT-convolution variant. Both show excellent agreement with theory. The sinc function for temporal coherence is somewhat poorly suited to the hybrid approach, as its tails retain significant value for large time differences. Here, only the values of the sinc function inside the first five side lobes are retained, making the filter (W) 219 weights (coefficients) long. Often, significantly fewer coefficients can be used. Even so, the run times and memory requirements of the standard and hybrid FFT-convolution approaches are quite different. Both simulations are executed using MATLAB® 2013A on a computer with dual Intel Xeon E5-2690 processors and 192 GBytes of random access memory (RAM). The standard implementation requires 181 s and 116 GBytes of RAM, while the hybrid FFT-convolution approach requires 9,160 s but only 0.7 GBytes of RAM. Of note, the hybrid's memory requirement will not increase with the total number of time samples, as it depends only on the spatial samples and filter coefficients. Thus, the hybrid FFTconvolution method greatly reduces memory requirements at the expense of computation time.
At this point, it is worth noting that the hybrid FFT-convolution method developed in this work offers tremendous practical benefits. The memory needs of the complex screen and phase screen methods can be very expensive. For example, Rydberg and Bengtsson use the complex screen method to measure the accuracy of an alternate numerical approach [40]. They do this for about 10 7 time samples but only a single point in space. If spatial variation is also important, such as for propagation through turbulence, then the field must be simulated in space as well as time. If 256 spatial points per side are needed, then 10 7 time samples will require about 38.7 TBytes of RAM using the standard complex screen method. Such computations are possible via distributed-memory FFT algorithms (e.g. FFTW) running on supercomputers. However, distributed-memory simulations require orders of magnitude more resources than runs on a conventional computer, and scheduling such simulations can be troublesome, as multiple research teams often compete for the same resources. On the other hand, the hybrid implementation would still require only about 0.7 GBytes of RAM, which is easily supplied by a conventional computer. Such a computer can readily be connected directly to laboratory hardware, as demonstrated in the next section. Notably, the time steps of the hybrid method are usually only limited by how long one is willing to wait for the simulation to complete. Thus, the hybrid method offers very significant practical benefits.

Laboratory young's experiment
In addition to the numerical simulations, a laboratory Young's experiment was performed using an SLM. The SLM is flood illuminated by a highly-coherent 532 nm laser, as shown in Fig. 4. Because the SLM can easily control phase, while control of both amplitude and phase is more difficult, the phase screen method is used to simulate the source. More specifically, the hybrid FFT-convolution variant is used, as the computer controlling the SLM cannot support sufficient memory for the other methods. The source phase is placed on the SLM at each simulated time step. Also, a grating phase is placed on the SLM such that all light striking the pinholes is directed into the desired first diffraction order. Finally, rather than using a physical lens, the lens phase is also added to the SLM.
The conditions of the experiment include Gaussian forms for the spatial and temporal coherence functions such that the DoC is given by Here, δ = 4.24 mm, while v Δ = 99.5 GHz. The pinholes are d = 1 mm in diameter with a center-to-center separation of D = 1.5 mm. Source sampling is defined by the 6.4 μm pixel pitch of the Holoeye LETO SLM. The experiment runs for a total of 500 coherence times with a time sampling of 0.667 ps (1/10 th the coherence time), and the field emitted from one of the pinholes is delayed by 1/2 the coherence time. The focal length is 2 m. The optical detector is a Lumenera Lw135RM with a 4.65 μm pixel pitch, which defines the observation plane sampling. Flat-fielding and dark frame subtraction are used to improve the quality of the results. Fig. 5 compares the experimental time-averaged irradiance to theory. The agreement is excellent with only minor differences in both the 2D irradiance images and the slice plots. Further, the associated video file, Visualization 1, shows the time-evolving irradiance for the first 100 coherence times of the experiment. It shows how the fringes shift over time due to phase changes across the pinholes. These changes are caused by the partial spatial and temporal coherence of the source. Observing the time evolution of the fringes provides physical insight into the formation of the time-averaged irradiance.
Fig. 5 appears to show very little noise. Potential sources of noise include pixel non-uniformity, dark current, stray light, read-out noise, and registration. Non-uniformity is greatly reduced by flat-fielding, while dark current is mitigated by dark frame subtraction. Spatial filtering of the first-order grating reflection and use of an enclosure help with stray light, while post-processing is used to estimate and remove the remainder. The post-processing simply involves measuring the mean irradiance at points five diffraction lobes away from the center of the fringe pattern. That mean value is then subtracted from the images, almost entirely removing the stray light. The read-out noise is largely removed by the time-averaging operation. Finally, registration effects are visible but small, indicated by the slight spatial shift between the experimental and theoretical fringe patterns in Fig. 5c.

Conclusion
In this paper, three methods were presented for the synthesis of continuous-wave or pulsed time-evolving fields with partial spatial and temporal coherence. The complex screen and phase screen approaches involve filtering in the frequency domain, allowing use of the numerically efficient FFT algorithm. The complex screen method permits a wide range of coherence functions, even coherence functions which are not separable in space and time, but its memory requirements are demanding. The phase screen method generates fields which are random in phase only, which is helpful for some applications, such as field generation using an SLM. However, producing partiallycoherent fields with non-exponential (i.e., not Gaussian or Lorentzian) DoCs can be difficult. The third method is a different implementation of the first two. It combines spectral filtering in space with convolution in time, greatly reducing the memory requirements at the expense of computation time. This implementation greatly eases synthesis of large sources by avoiding the need for distributed-memory supercomputers.
Simulations and experiments were performed to test the validity of all three synthesis methods. Using a modified Young's interference experiment, all methods showed excellent agreement with theory in both the numerical and laboratory experiments. These results demonstrated the accuracy of the three simulation methods.
The work presented here will be useful in applications where the time evolution of a partially-coherent source is important. These applications include scenarios involving high-speed electronics/detec- Fig. 4. Diagram of the Young's experiment. The SLM is flood illuminated by a highly coherent 532 nm laser. A grating phase is applied to the SLM region representing the pinholes. The first-order reflection passes through a 4-f optical system with an iris in the focal plane to spatially filter out the undesired diffraction orders. The field then propagates to the detector. tors, propagation of optical vortices through turbulence, et cetera. In addition to its practical uses, this work can also be used as an educational tool, in which watching the source temporally evolve can provide more insight than theory alone.