A numerical wave optics approach for simulating a partial spatially coherent beam is presented. The approach involves the application of a sequence of random phase screens to an initial beam field and the summation of the intensity results after propagation. The relationship between the screen parameters and the spatial coherence function for the beam is developed and the approach is verified by comparing results with analytic formulations for a Gaussian Schell-model beam. The approach can be used for modeling applications such as free space optical laser links that utilize partially coherent beams.
©2006 Optical Society of America
The partially (spatial) coherent beam (PCB) has been studied extensively over the last four decades as an approach for improving the performance of various laser system applications including free-space optical laser communications [1,2,3]. The characteristics of PCBs and methods for producing such beams have been examined in a number of publications [4,5,6]. Most studies consider the Gaussian Schell-model (GSM) beam , which is an analytically tractable model where the beam field amplitude distribution and the spatial coherence function are both Gaussian. However, practical PCBs may not follow the GSM theory and are better examined in detail through some type of numerical simulation approach such as a wave optics simulation. An example is the pseudo-partially coherent beam (PPCB) suggested for use in communications applications [7,8]. In addition, explicit modeling of the optical systems and propagation scenarios associated with PCB applications may be best done through wave optics simulations. For these reasons we have developed an approach for modeling the PCB in a wave optics simulation. The approach can accommodate arbitrary initial beam amplitude and spatial coherence functions.
In this paper we describe the modeling approach and relate the method to the Schell theorem that describes the Fraunhofer intensity pattern of a PCB. To validate the simulation approach, we compare intensity patterns produced by the method with analytical results for a nominal GSM beam propagation and a GSM beam incident on a two-pinhole aperture.
2. Simulation approach
The simulation approach is outlined as follows: a spatially random phase screen is applied to a deterministic beam field, the result is propagated and the intensity is formed. The process is repeated many times with different realizations of the phase screen and the resulting intensity patterns are averaged to produce the PCB intensity. Analytically, we can consider the sequence of screens to approximate a phase function that is spatially random and time-varying and the received intensity as a time-averaged quantity. The method is analogous to an approach for modeling the effects of temporal coherence involving a spectral intensity average .
Consider the field described by
where U 0(p, q) is the deterministic field of an initial beam defined at spatial position (p, q) and tA is a transmittance function given as
where ξ (p, q; t) is a random phase introduced at time t. After propagation, the time-averaged intensity of the Fraunhofer pattern is
where z is the propagation distance, λ is the wavelength and the angular brackets with subscript t denote a time average. Equation (3) is representative of the process implemented in our simulation approach. Expanding the squared-modulus quantity and interchanging the averaging and integration operations we have
We introduce the following average and difference notation p¯ = (p 1 + p 2)/2, q¯ = (q 1 + q 2)/2, Δp = p 1 - p 2, Δq = q 1 - q 2. Assuming the random phase is stationary in the first increment, the correlation of the transmittance function R(p 1,q 1,p 2,q 2;t) = 〈exp(jξ(p 1,q 1;t))exp(jξ(p 2,q 2;t))〉, depends only on the coordinate differences Δp, Δq. As a result, R(p 1, q 1, p 2, q 2; t) = R(Δp, Δq; t) and Eq. (4) can be expressed as 
where U 0 is the autocorrelation function of the deterministic beam U 0
Assuming the time-average is equivalent to an ensemble average, in Eq. (5), we could write 〈R(Δp,Δq;t)〉t = R(Δp,Δq). Equation (5) is a statement of Schell’s theorem , which describes the Fraunhofer intensity profile for a partial spatially coherent beam with a spatial coherence function corresponding in this case to R(Δp,Δq). Thus our modeling approach is consistent with Schell’s theorem.
In the numerical simulation, the time varying transmittance function is approximated by a sequence of independent, random phase screens tA (p,q)i =exp(jξ(p,q)i) where i is an integer i∊ [1,K]. For each realization of the transmittance function, the propagation of the deterministic field-transmittance product U 0(p,q)∙tA (p,q) can be accomplished using a standard numerical propagation technique such as an FFT-based approach. A K-average of the resulting intensity frames creates the simulated intensity profile of the PCB.
We note that the transmittance function in Eq. (2) could also be an absorbing screen that affects the magnitude of the deterministic beam. With the appropriate function, this type of screen can also be used to simulate a partially coherent beam. However, the random phase screen has the advantage of preserving the total power of the deterministic beam.
3. Phase screen for Gaussian Schell-model beam
The key to applying the simulation technique to a PCB of interest is defining the phase screen transmittance function since the transmittance autocorrelation function is equivalent to the beam’s coherence function. For illustration purposes we discuss the GSM beam, which has a Gaussian coherence function . We begin by modeling the random phase function as
where ⊗ indicates a spatial convolution, r(q, p; t) is a spatially uncorrelated or “delta correlated” random signal with standard deviation σr and f(p, q) is a Gaussian function that acts as a filter with normalized area and standard deviation σf
The parameter σr is related to the amplitude variation of the screen, something like optical path length variation, and σf is a transverse spatial correlation length parameter. Both parameters have units of length. The convolution of Eq. (7) results in a spatially correlated random function. Following the approach by Goodman , we obtain the transmittance correlation function
where Δξ = ξ(p 2,q 2;t)-ξ(p 1,q 1;t) and note that a multiplicative phase term exp(j〈Δξ〉t) has been dropped from the right side of Eq. (9) because 〈Δξ〉t = 〈ξ(p 2,q 2;t)〉t - 〈ξ(P 1,q 1;t)〉t = 0. is the variance of Δξ, and
where Γξ(∙) indicates the spatial autocorrelation function of the phase ξ, and can be taken to be 
Again, we write Rrr (Δp, Δq) = 〈Rrr (Δp, Δq;t)〉 as the autocorrelation function of r(q, p; t), which is found to be
where δ(Δp, Δq) is the delta function. The second part of the right side of Eq. (11) is derived as
Consequently, the transmittance autocorrelation function is
We note if /(4π∙) ≫ 1, then Eq. (15) can be approximated as
where the variance of this Gaussian function is
Thus, Eq. (16) provides a Gaussian correlation function as required for the GSM beam.
For implementation in a wave optics simulation, the GSM transmittance phase screen can be created by convolving two arrays defined by the discrete functions
where n and m are integers, N and M are the dimensions of the array (number of samples) in the horizontal and vertical directions, and, Δn = Ln /N and Δm = Lm /M represent the physical spacings between the samples where Ln and Lm are the physical lengths of the horizontal and vertical sides of the array. rand(n,m) denotes a random number function with variance σr that provides a value for each (n, m) position. The division by (ΔnΔm)1/2 in Eq. (19) is required for conversion to the unit-less sample domain. The convolution is easily accomplished in the computer using the Fourier convolution theorem by Fourier transforming the functions in Eqs. (18) and (19), multiplying the results point-wise, inverse transforming and finally multiplying by Δn∙Δm to correctly approximate the convolution integral.
Note there is some flexibility in the choice of σr and σf when modeling the GSM beam. A range of combinations can yield a specific coherence size parameter [Eq. (17)] and satisfy the approximation required in going from Eq. (15) to Eq. (16). We found in practice the choice of a relatively larger value of σf tends to produce a smoother intensity pattern with fewer number of frames. However, aliasing can become a problem when σf is chosen too large where the filter function extends significantly beyond the effective boundaries of the numerical array.
4. Simulation examples
To verify the wave optics approach, we compared simulation results with analytic results for a GSM beam. The analytic expressions for the GSM beam are presented here for reference. For a Gaussian transmitted beam field with initial peak value of 1 and radius w 0 (1/e field value), after propagating a distance z, the optical intensity (irradiance) can be expressed as 
where k = 2π/λ is the wave number, λ is the nominal wavelength, and is the variance of the Gaussian coherence function, which corresponds to the parameter defined in Eq. (17). Figure 1 shows a comparison of simulation and analytic intensity patterns for two sets of parameters. The array size was 256×256 points and a uniform random number generator was used for the screens. The propagation operation in the simulation was implemented in AOTools, an adaptive optics software toolbox for MatLab. 2000 phase screen realizations (K= 2000) were used for the simulation results. The simulation and analytic profiles are identical except for small random variations in the simulation intensity due to the finite number of frames. The simulation profiles continue to become smoother as more frames are averaged.
In a second test, we simulated the intensity pattern of a GSM beam that is initially incident on a mask of two small holes. For holes of diameter d that are equidistance from the optical axis and whose separation is much smaller than the incident beam radius, the normalized Fraunhofer intensity is approximately given by 
where J 1(.) is the first-order Bessel function of the first kind and μ is the complex coherence factor where
and Δṕ and σq́ are the pinhole spacings in the horizontal and vertical directions, respectively. Figure 2 shows a comparison of simulation and analytic fringe profiles for a Gaussian beam of radius 8 cm incident on two pinholes of diameter of 3 mm and spacing of 1 cm. μ ranges from 0.15 to 0.6 for the different cases. 200 realizations with array size 512×512 points were used and the simulation and analytic curves again agree nearly completely except for a slight contrast difference in the first case and variations in the simulation results due to the finite number of realizations.
A multi-phase screen approach has been described for modeling the propagation of a PCB in a wave optics simulation. The approach was shown to be consistent with Schell’s theorem and results were compared with analytic theory for several cases involving Gaussian Schell-model beams. The comparisons showed the simulation intensity patterns are nearly identical to the analytic patterns except for variations due to finite frame averaging.
An important consideration for the approach is the number of frames used to generate a PCB intensity pattern. The pattern approaches the theoretical result with more frames but the simulation run time can become prohibitive. As illustrated in the second example, several hundred realizations may be adequate for some GSM beam cases. However, in future work, a metric such as mean squared error could be used to study how the numerical results converge to the theory as realization number K increases. Other types of phase screens, for example, quasi-random screens, might be explored to improve the behavior of the metric and result in faster execution times. These investigations would help to advance the practical implementation of the method.
This manuscript was prepared through collaborative participation in the Communications and Networks Consortium sponsored by the U. S. Army Research Laboratory under the Collaborative Technology Alliance Program, Cooperative Agreement DAAD19-01-2-0011.
References and links
1. S. C. H. Wang and M. A. Plonus, “Optical beam propagation for a partially coherent source in the turbulent atmosphere,” J. Opt. Soc. Am. 69, 1297–1304 (1979). [CrossRef]
2. Y. Baykal, M. A. Plonus, and S. J. Wang, “The scintillations for a weak atmospheric turbulence using a partially coherent beam,” Radio Sci. 18, 551–556 (1983). [CrossRef]
3. J. C. Ricklin and F. M. Davidson, “Atmospheric turbulence effects on a partially coherent Gaussian beam: implication for free-space laser communication,” J. Opt. Soc. Am. A 19, 1794–1802 (2002). [CrossRef]
4. W. Martienssen and E. Spiller, “Coherence and fluctuations in light beams,” Am. J. Phys. 32, 919–926 (1964). [CrossRef]
5. H. Arsenault and S. Lowenthal, “Partial coherence of an object illuminated with laser light through a moving diffuser,” Opt. Commun. 1, 451–453 (1970). [CrossRef]
6. E. Tervonen, A. T. Friberg, and J. Turunen, “Gaussian Schell-model beams generated with synthetic acousto-optic holograms,” J. Opt. Soc. Am. A 9, 796–803 (1992). [CrossRef]
7. D. G. Voelz and K. J. Fitzhenry, “Pseudo-partially coherent beam for free-space laser communication,” in Free-Space Laser Communications IV, J. C. Ricklin and D. G. Voelz, eds., Proc. SPIE 5550, 218–224 (2004). [CrossRef]
8. X. Xiao and D. G. Voelz, “Wave optics simulation of partially coherent beams,” in Free-Space Laser Communications V, D. G. Voelz and J. C. Ricklin, eds., Proc. SPIE 5892, 219–227 (2005).
10. J. W. Goodman, Statistical Optics, (John Wiley & Sons, 1985).
11. L. Mandel and E. Wolf, “Radiation from sources of any state of coherence,” in Optical Coherence and Quantum Optics, (Cambridge University, 1995), pp. 229–337.
12. H. Stark and J. W. Woods, Probability and Random Process with Applications to Signal Processing, (Prentice Hall, 2002), Chap. 7.