Numerical Approach for Studying the Evolution of the Degrees of Coherence of Partially Coherent Beams Propagation through an ABCD Optical System

Abstract

In this paper, we propose a numerical approach to simulate the degree of coherence (DOC) of a partially coherent beam (PCB) with a Schell-model correlator in any transverse plane during propagation. The approach is applicable for PCBs whose initial intensity distribution and DOC distribution are non-Gaussian functions, even for beams for which it is impossible to obtain an analytical expression for the cross-spectral density (CSD) function. Based on our approach, numerical examples for the distribution of the DOC of two types of PCBs are presented. One type is the partially coherent Hermite–Gaussian beam. The simulation results of the DOC agree well with those calculated from the analytical formula. The other type of PCB is the one for which it is impossible to obtain an analytical expression of CSD. The evolution of the DOC with the propagation distance and in the far field is studied in detail. Our numerical approach may find potential applications in optical encryption and information transfer.

1. Introduction

Over the past decades, partially coherent beams (PCBs) have attracted wide attention due to their important applications, such as particle trapping/manipulation, wide-field optical coherence tomography, laser material processing, optical metrology, 3D imaging and, so on [1,2,3,4,5,6,7,8,9]. Generally, PCBs are characterized by the so-called cross-spectral density (CSD) function in the space–frequency domain or mutual coherence function (MCF) in the space–time domain, and the coherence characteristics of PCBs are described by the degree of coherence (DOC) which can be obtained from the CSD function or the MCF. It is known that the DOC plays a central role in determining the beam propagation characteristics, image resolution, and effects of interaction with matter [10,11,12,13,14]. However, most studies are focused on the PCBs whose DOCs are of Gaussian function because it is a hard task to devise a physically realizable CSD function with arbitrary DOC. Recently, Gori and co-workers proposed a sufficient condition for mathematically devising a genuine CSD function [15,16]. Since then, researchers have turned their attention to study PCBs with special DOCs (non-Gaussian function) [17,18,19,20,21]. Owing to their special DOCs, these beams exhibit peculiar propagation features even in free space. For example, a non-uniformly correlated (NUC) beam displays “self-focusing” and lateral shift of intensity maxima during free-space propagation [17]. The multi-Gaussian-correlated beam and Laguerre–Gaussian (LG)-correlated beam will produce a flat-topped beam profile and ring-shaped beam profile in the far field, respectively [18,19]. In addition, the DOCs may carry the information of the beam itself or of the object with which the PCBs interact [22]. For example, by measuring the far-field DOC [23], the azimuthal and radial mode orders of a partially coherent LG_pl beam can be determined. In [14], the authors point out that, when a PCB is obstructed by an opaque obstacle, its DOC can self-reconstruct. In addition, several methods have been proposed to experimentally generate the prescribed DOC of a PCB, for example by using a spatial light modulator, intra-cavity modulation and a deformable mirror, etc. [24,25,26,27]. However, the aforementioned studies are confined to PCBs whose initial intensity distributions have a Gaussian function.

On the other hand, PCBs with both initial non-Gaussian beam profiles and DOCs have also received considerable attention in recent years [28,29,30,31,32,33,34]. In the following text, we use the term general PCBs (GPCBs) to include a class of PCBs with non-Gaussian amplitudes and DOCs. Through controlling the DOCs and beam profiles in the source plane simultaneously, GPCBs can generate a far-field beam array with controllable forms [32,33], and the profile of each beamlet in the array is also adjustable, making them ideal sources for particle trapping and manipulation. Nevertheless, in these studies, the evolution properties of the DOC during propagation were neglected. Though other studies [35,36] presented the reciprocity relation between the initial beam profile and the DOC in the far filed, this is only applicable for special kinds of GPCBs, known as quasi-homogeneous beams. In such beams, the variation of the DOCs with the spatial point r is much faster than that of intensity distribution with r. Up to now, only few papers have dealt with the evolution of the DOCs during GPCBs passage through free space or paraxial optical systems, especially for those beams whose DOCs or initial intensity distributions are impossible to obtain using an analytical expression.

In this paper, we propose a numerical approach for calculating the DOC of a GPCB in any transverse plane perpendicular to the propagation axis during propagation in free space or paraxial optical systems. In this approach, we introduce white noises [37,38,39] into the convolution approach [40] to generate speckle patterns in any transverse plane of propagation. Then, the DOC is retrieved from the fourth-order correlation function of the speckle patterns. On the basis of this approach, we investigate the effects of source intensity distribution and coherence width on the DOC during propagation in free space through several numerical examples.

2. Theory

Let us consider a scalar, quasi-monochromatic, and statistically stationary optical beam, propagating along the z-direction. In the space–frequency domain, the second-order statistics of such GPCB are typically characterized by the CSD function $W({\mathsf{\rho }}_1,{\mathsf{\rho }}_2,z)=\langle E^{*}({\mathsf{\rho }}_1,z)E({\mathsf{\rho }}_2,z)\rangle$, where E denotes the random electric field, ${\mathsf{\rho }}_1$ and ${\mathsf{\rho }}_2$ are two arbitrary position vectors in the transverse plane, the asterisk is the complex conjugate, and the angular brackets is the ensemble average over the field fluctuations. Assume that the GPCB we consider is a Schell-model beam whose DOC function is only dependent on the separation of two points. In this situation, the CSD function in the source plane (z = 0) can be written as [15,16]

$$W_0\left({\mathbf{r}}_1,{\mathbf{r}}_2\right)={\tau }_0\left({\mathbf{r}}_1\right){\tau }_0^{*}\left({\mathbf{r}}_2\right){\mu }_0\left({\mathbf{r}}_1-{\mathbf{r}}_2\right),$$

(1)

where ${\mathbf{r}}_1$ and ${\mathbf{r}}_2$ are two arbitrary position vectors in the source plane, ${\tau }_0$ is a complex amplitude function, and ${\mu }_0$ is the DOC function. For a physically realizable CSD function, the DOC can be alternatively represented the following integral form [19]

$${\mu }_0\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)={\displaystyle \int W_s\left({\mathbf{v}}_1,{\mathbf{v}}_2\right)}\exp \left[-i2\pi \left({\mathbf{r}}_1-{\mathbf{r}}_2\right)\cdot \mathbf{v}\right]d^2\mathbf{v},$$

(2)

with

$$W_s\left({\mathbf{v}}_1,{\mathbf{v}}_2\right)=\sqrt{p\left({\mathbf{v}}_1\right)}\sqrt{p\left({\mathbf{v}}_2\right)}\delta \left({\mathbf{v}}_1-{\mathbf{v}}_2\right).$$

(3)

The function p(v) in Equation (3) must be non-negative for any v to ensure that the CSD is genuine; $\delta$ denotes the Dirac function. In order to introduce our numerical approach, we rewrite the $\delta$ function in the follow form

$$\delta \left({\mathbf{v}}_1-{\mathbf{v}}_2\right)=\langle C_n\left({\mathbf{v}}_1\right)C_n^{*}\left({\mathbf{v}}_2\right)\rangle ,$$

(4)

where C_n(v) (n = 1, 2, 3, …) denotes the random complex function whose probability density functions of amplitude and phase of C_n obey the negative exponent and uniform distribution, respectively [39]. The angular bracket denotes the ensemble average. Here, C_n(v) can be viewed as the field of the white noise.

By substituting Equation (4) into Equations (1)–(3) and after some mathematical operations, Equation (1) turns out to be

$$W_0\left({\mathbf{r}}_1,{\mathbf{r}}_2\right)=\langle T_{n0}\left({\mathbf{r}}_1\right)T_{n0}^{*}\left({\mathbf{r}}_2\right)\rangle \approx \frac{1}{N}{\displaystyle \sum _{n0=1}^N{T}_{n0}\left({\mathbf{r}}_1\right)T_{n0}^{*}\left({\mathbf{r}}_2\right)},$$

(5)

with

$$T_{n0}\left({\mathbf{r}}_i\right)={\tau }_0\left({\mathbf{r}}_i\right)\times F_T\left[\sqrt{p\left({\mathbf{r}}_i\right)}{C}_n\left({\mathbf{r}}_i\right)\right],i=1,2.$$

(6)

where F_T stands for the symbol of Fourier transform, T_n0 (r) (n = 1, 2, 3…) denotes the instantaneous electric fields, each of which represents one realization of a GPCB. In Equation (5), we use the incoherent superposition of instantaneous random realizations instead of ensemble averaging over the random fields, since we assume the beam is statistically stationary. If the number of realizations N is large enough, the incoherent realization superposition is a good approximation. The first term ${\tau }_0$ in Equation (6) is the complex amplitude, as we described above, and the second term contains the information of the DOC of the source. To synthesize the realization of T_n0 (r) numerically, one can use commercial software such as MATLAB to generate the random complex function C_n and then take the Fourier transform of the product of the function C_n and the square root of the p(v) function. Finally, one instantaneous electric filed is acquired by multiplying the Fourier result by the function ${\tau }_0$. Note that the DOC is determined by the function p(v) in the simulation. Once each instantaneous electric field in the source plane is known, the relation between such fields in the source plane and in the output plane passing through a paraxial ABCD optical system is linked to the Collins integral formula, given by

$$T_{nz}\left(\mathsf{\rho }\right)=\left(-\frac{i}{\lambda B}\right){\displaystyle \int T_{n0}\left(\mathbf{r}\right)}\exp \left[\frac{ik}{2B}\left(A{\mathbf{r}}^2-2\mathsf{\rho }\cdot \mathbf{r}+D{\mathsf{\rho }}^2\right)\right]d^2\mathbf{r}，$$

(7)

where $k=2\pi /\lambda$ is the wavenumber associated with $\lambda$, which is the wavelength of the light, and $\mathsf{\rho }$ denotes the arbitrary position vector in the output plane. A, B, D are the elements of the transfer matrix of an optical system. On Substituting Equation (6) into (7), we obtain the expression

$$T_{nz}\left(\mathsf{\rho }\right)=\left(-\frac{i}{\lambda B}\right)\exp \left(\frac{ikD{\mathsf{\rho }}^2}{2B}\right){\displaystyle \int F_T\left[\sqrt{p\left(\mathbf{r}\right)}{C}_n\left(\mathbf{r}\right)\right]{\tau }_0\left(\mathbf{r}\right)\exp \left(\frac{ikA{\mathbf{r}}^2}{2B}\right)\times }\exp \left(\frac{ik\mathsf{\rho }·\mathbf{r}}{B}\right)d^2\mathbf{r}.$$

(8)

After some rearrangement, Equation (8) can be written in the simple convolution form

$$T_{nz}\left(\mathsf{\rho }\right)=\left(-\frac{i}{\lambda B}\right)\exp \left(\frac{ikD{\mathsf{\rho }}^2}{2B}\right)\left(\sqrt{p}{C}_n\right)\otimes \tilde{U}\left(-\frac{\mathsf{\rho }}{\lambda B}\right)，$$

(9)

where $\tilde{U}$ is the inverse Fourier transform of the new defined function $U(\mathbf{r})={\tau }_0\left(\mathbf{r}\right)\exp \left(ikA{\mathbf{r}}^2/2B\right)$, and $\otimes$ denotes the convolution symbol. Equation (9) establishes the relation of an instantaneous field between the source plane and the output plane with the convolution form, which is the main result of this paper. It provides one convenient way to simulate the second-order statistics of a GPCB in the output plane.

Here, we describe the method for numerical simulating the instantaneous electric field $T_{nz}\left(\mathsf{\rho }\right)$ using MATLAB. It follows from Equation (9) that the key point is to calculate the convolution term $\left(\sqrt{p}{C}_n\right)\otimes \tilde{U}$. According to the convolution theorem, such term can be expressed as $F_T^{-1}\left[U\times F_T^{}\left(\sqrt{p}{C}_n\right)\right]$, where $F_T^{-1}$ denotes the symbol for inverse Fourier transform. The operations for the $F_T^{-1}$ and $F_T$ in MATLAB are based on the fast Fourier transform (FFT) algorithm. In numerical simulation, each realization of any function is represented as a matrix with N₁ × N₁ number of elements. We first assume the function U has the physical dimensions L × L, where L is the length of the coordinate. The corresponding coordinate r = (x, y) for x or y is from −L/2 to L/2 with step L/N₁. In order to match the coordinate (x, y), the function p (related to the DOC) must sample in the spatial frequency domain having dimensions N₁/L × N₁/L, and the corresponding coordinate f_x or f_y is from −N₁/2L to N₁/2L, with step 1/L [41]. The random complex function C_n is generated by [randn(N₁) + i × randn(N₁)]/sqrt(2), where randn(N₁) is the build-in function returning the N₁ × N₁ matrix consisting of random numbers of normal distribution. Then, the Fourier transform function $\sqrt{p}{C}_n$ is multiplied by the function U, and, finally, $T_{nz}\left(\mathsf{\rho }\right)$ is obtained by multiplying the inverse transform of the above results $U\times F_T^{}\left(\sqrt{p}{C}_n\right)$ by the coefficient before the convolution shown in Equation (9). The corresponding coordinates for ρ_x or ρ_y in the out plane will be from −λBN₁/2L to λBN₁/2L, with step λB/L [41].

According to Equation (9), the average spectral density is calculated by the following formula

$$S_z\left(\mathsf{\rho }\right)={\displaystyle \sum _{n=1}^N{S}_{nz}\left(\mathsf{\rho }\right)}/N,$$

(10)

where $S_{nz}(\mathsf{\rho })=T_{nz}^{*}(\mathsf{\rho })T_{nz}(\mathsf{\rho })$ is the instantaneous spectral density (one realization).

In order to evaluate the DOC, we first calculate the normalized intensity correlation (fourth-order correlation function) between points ${\mathsf{\rho }}_1$ and ${\mathsf{\rho }}_2$, i.e.,

$$g^{\left(2\right)}\left({\mathsf{\rho }}_1,{\mathsf{\rho }}_2\right)=\frac{{\displaystyle \sum _{n=1}^N{S}_{nz}({\mathsf{\rho }}_1)S_{nz}({\mathsf{\rho }}_2)}}{N{S}_z({\mathsf{\rho }}_1)S_z({\mathsf{\rho }}_2)}.$$

(11)

By applying the Gaussian moment theorem [1], the modulus of square of the DOC can be extracted from the normalized intensity correlation function, given by

$${\left|{\mu }_z\left({\mathsf{\rho }}_1,{\mathsf{\rho }}_2\right)\right|}^2=g^{\left(2\right)}\left({\mathsf{\rho }}_1,{\mathsf{\rho }}_2\right)-1.$$

(12)

It follows from Equation (12) that one can numerically calculate the DOC in the output plane from the normalized intensity correlation function. Note that Equation (12) is only valid for the GPCBs with a Schell-model type. Hence, the DOC is only a function of the separation of two points ${\mathsf{\rho }}_1$ and ${\mathsf{\rho }}_2$, i.e., ${\mu }_z\left({\mathsf{\rho }}_1,{\mathsf{\rho }}_2\right)={\mu }_z\left({\mathsf{\rho }}_2-{\mathsf{\rho }}_1\right)$.

3. Simulation Results

In this section, to demonstrate the validity for the numerical approach, we first compare the results for the DOC obtained with our numerical approach with those derived from an analytical expression. Then, we present numerical examples for the evolution of the DOC with the propagation distance, using our numerical method. In this numerical case, it is impossible to obtain the analytical expression of the GPCB during propagation. In the following simulation, the number of elements N₁ × N₁ and N are chosen to be 512 × 512 and 10⁴, respectively. Under this circumstance, it takes about 3.7 min for each plot. The computer is from Dell, Texas, USA, whose configurations used in the simulation are i7-8700 CPU @ 3.2GHz and 16GB RAM. Here, we want to stress that the advantage of our numerical approach is that it can deal with partially coherent beams that cannot be described by an analytical expression.

In order to demonstrate the validity of our simulation method, a typical Schell-model beam, named a partially coherent elegant Hermite–Gaussian (PCEHG) beam [42], is used to compare the evolution of DOC from the analytical results with that from the simulation results. The CSD function of the PCEHG beam in the source plane is expressed as

$$\begin{array}{l}{W}_0\left({\mathbf{r}}_1,{\mathbf{r}}_2\right)={\tau }_0\left({\mathbf{r}}_1\right){\tau }_0^{*}\left({\mathbf{r}}_2\right){\mu }_0\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)\\ =H_n\left(\frac{x_1}{{\omega }_0}\right)H_m\left(\frac{y_1}{{\omega }_0}\right)H_n\left(\frac{x_2}{{\omega }_0}\right)H_m\left(\frac{y_2}{{\omega }_0}\right)\\ \times \exp \left(-\frac{{\mathbf{r}}_1^2+{\mathbf{r}}_2^2}{{\omega }_0^2}\right)\exp \left(-\frac{{\left|{\mathbf{r}}_1-{\mathbf{r}}_2\right|}^2}{2{\sigma }_g^2}\right).\end{array}$$

(13)

$H_n$ is the Hermite polynomial of order n, ${\omega }_0$ and ${\sigma }_g$ are the beam waist size and the transverse coherence width, respectively. The analytical expression for the CSD function of the PCEHG beam at distance z in free-space propagation can be found in [42] (see Equation (17) in [42], where the elements of the transfer matrix are A = 1, B = z, and D = 1). Figure 1a1–e1 illustrates the evolution of the modulus of the square of the DOC with the propagation distance in free space calculated from the analytical formula. The beam parameters in the calculation are chosen to be n = 2, m = 0, $\lambda =532\mathrm{nm}$, ${\omega }_0=1\mathrm{mm}$, and ${\sigma }_g=1\mathrm{mm}$. For comparison, the corresponding results attained by our numerical approach are shown in Figure 1a2–e2. It can be seen that the simulation results agree reasonably well with those obtained from the analytical formula, implying that our method is valid and effective in simulating the evolution of the DOC during propagation.

As we state in Section 2, the advantage of the numerical approach is that it can deal with the evolution of the DOC of a GPCB whose CSD function during propagation is impossible to acquire by an analytical expression. In the following example, we assume that the function p(v) (Fourier transform of the DOC) is a soccer pattern, as shown in Figure 2a; the corresponding modulus of the square of the DOC is plotted in Figure 2c. The initial amplitude function ${\tau }_0$ we choose is the elegant Hermite–Gaussian distribution, as follows

$${\tau }_0\left(x,y\right)=H_n\left(x/{\omega }_0\right)H_m\left(y/{\omega }_0\right)\exp \left[-\left(x^2+y^2\right)/{\omega }_0^2\right].$$

(14)

Figure 2b,d presents the source intensity distribution of such GPCB and its modulus of the square of the Fourier transform, respectively. The relevant parameters are set as n = 2, m = 0, ${\omega }_0=1\mathrm{mm}$, and $\lambda =532\mathrm{nm}$. If the source beam is assumed to be a quasi-homogeneous planar source [1,35,36], then, Figure 2d can be treated as the modulus of the square of the DOC in the far field.

Before presenting the simulation results, we want to determine the relationship between the initial coherence width and the size of the soccer pattern. As stated in [1,10,43], the source transverse coherence width is inversely proportional to the size of the function of p. Here, we assume that the transverse coherence width is ${\delta }_0=1/l$ where l denotes the size of the distribution of p(v) (radius of the soccer).

Figure 3 presents the simulation results of the DOC of the GPCB with different initial coherence width at several propagation distances. The parameters used in the simulation are n = 2, m = 0, ω₀ = 1 mm, and λ = 532 nm. One can see that the DOC pattern can keep the profile similar to that in the source plane (see in Figure 2c), in a range from z = 0 to certain propagation distances, but this range is closely dependent on the initial coherence width. When the coherence width is large (${\delta }_0=0.64\mathrm{nm}$), the DOC can keep its profile invariant at about 0.7 m (see Figure 3a3–c3]. In the case of ${\delta }_0=0.32\mathrm{nm}$, the profile of the DOC changes much faster compared to that in the case for ${\delta }_0=0.64\mathrm{nm}$ and evolves into a shape similar to that in the far field at distance z = 1.0 m (see Figure 2d and Figure 3d1).

Let us now turn our attention on the influence of the initial coherence width on the DOC pattern of such GPCB in the far field. Figure 4a–c shows our numerical simulations of the DOC distributions at propagation distance z = 1000 m (far field) with different values of coherence width. The simulation parameters are the same as those used in Figure 3. For the case of ${\delta }_0=0.32\mathrm{nm}$ (see Figure 4a), the DOC pattern is similar to that shown in Figure 2d, implying that the beam can be be considered as a quasi-homogeneous beam whose DOC in the far field is the Fourier transform of the source intensity distribution, independent of the DOC in the source plane. As the initial coherence width increases, the effects of the source DOC gradually appear. In this situation, the DOC pattern is determined by both the intensity distribution and the DOC function in the source plane. As the coherence width increases, the area of the side lobes becomes large. The enlargement of the side lobes in Figure 4b,c can be explained by the fact that, as the coherence width increases, the globe coherence in our calculation area, which is 2.0 mm × 2.0 mm, increases. As a result, the side lobes are enlarged in the high-coherence case (see in Figure 4b,c) compared with the low-coherence case (see in Figure 4a). It is also found that the distributions of the DOC in Figure 4 are left–right asymmetric, which may be caused by the fact that the p function (related to the DOC) has a slight left–right asymmetry, as shown in Figure 2a. To confirm our results, we choose as a circular function a p function, which is exactly circularly symmetric. The expression for a circular function is given by

$$p(\mathbf{v})=\{\begin{array}{l}1,\left|\mathbf{v}\right|\le l\\ 0,\left|\mathbf{v}\right|>l\end{array}$$

(15)

where l is the radius of soccer defined previously. Figure 4d–f presents the corresponding numerical simulations of the DOC distributions at distance z = 1000 with the corresponding coherence width. The p function is shown in Equation (15). One can see that the intensities of the side lobes in the left side and the right side are the same.

4. Conclusions

In summary, we have introduced a numerical approach to calculate the DOC distribution of a Schell-model type GPCB in any transverse plane perpendicular to the propagation axis during free-space propagation. This method first represents the GPCB as a sequence of random electric fields by introducing white noises into it. In this representation, each realization of random fields obeys the propagation law described by Huygens–Fresnel integral formula. Then, we use a convolution approach to obtain the corresponding realization of random fields in the output plane. At last, the information on DOC is acquired through evaluating the fourth-order correlation function of random fields in the output plane. The prominent advantage of this numerical approach is that it can conveniently calculate the DOC of any partially coherent Schell-mode beam, even when it is impossible to obtain using analytical expressions for the CSD function. Two numerical examples of the evolution of DOCs for two types of GPCBs with propagation distance were investigated in detail. One demonstrated that our approach is valid, by comparing our results with those obtained with the analytical formula. The other example showed the influence of the initial field amplitude and correlation function on the evolution of the DOC. The analytical formula for the CSD function cannot be applied to the beam in the second example. Our method can also deal with the DOC of a GPCB propagation through a paraxial ABCD optical system or atmospheric turbulence. We believe that this numerical approach may find applications in optical communication and optical encryption.