Electromagnetic Hanbury Brown and Twiss Effect in Atmospheric Turbulence

: The evolution of the 4 × 4 matrix with elements being the scintillation indices of the single-point Stokes parameters of a stationary electromagnetic beam-like optical ﬁeld in classic, weak atmospheric turbulence is revealed. It is shown that depending on the choice of the source parameters, the source-induced changes in the matrix elements of the propagating beam and those produced by turbulence can be either range-separated or conjoined. For theoretical analysis, the uniﬁed theory of coherence and polarization is used together with the extended Huygens-Fresnel integral approach. The results can be of interest for building robust communication and sensing systems operating in the presence of atmospheric ﬂuctuations.


Introduction
The discovery by Hanbury Brown and Twiss of non-trivial correlations in the intensity of spatially separated, fluctuating, classic light fields [1,2] has led to a number of groundbreaking technologies, the most important of which are the stellar interferometry [3] and ghost imaging [4]. In the majority of scenarios, the measured light field is thermal, i.e., obeys the circular Gaussian statistics, and, hence, involves a very simple relation between the field correlations and the intensity correlations, as a consequence of the Gaussian moment theorem [5].
Another area greatly benefiting from the optical intensity correlation analysis is the laser beam interaction with atmospheric turbulence, including a variety of applications, such as remote sensing, classic imaging, and free-space optical communications [6]. The laser beams, initially fully spatially coherent, are modulated by the air's fluctuating refractive index, being a very complex and non-Gaussian random process. The intensity correlations of the laser field interacting with turbulence hence carry a complex signature of the atmospheric refractive index's random process, which cannot be directly related to that acquired by field correlations.
Light fields radiated by thermal and partially coherent sources have also been studied on interaction with atmospheric turbulence, mostly under scalar treatment [7][8][9][10][11][12]. In such situations, the fluctuations in the optical field propagating through turbulence are also governed by a random process different from Gaussian. The usefulness of the spatially pre-randomized optical fields, however, is in the reduction of their intensity fluctuations as compared with those in a beam with the same geometry launched into the same turbulence channel by a laser. This idea was employed a number of times for improving the quality of free-space optical communications [13][14][15].
The original Hanbury Brown and Twiss experiment was based on the unpolarized, thermal light and relied on the scalar coherence theory. However, partially coherent light sources were later introduced that can have arbitrary polarization distribution and arbitrary vectorial spatial correlations [16,17]. Such sources can exhibit changes in coherence

The Average Stokes Parameters in Turbulence
The cross-spectral density matrix of an electromagnetic beam at positions r 1 and r 2 (where r 1 = (ρ 1 , z) and r 2 = (ρ 2 , z), ρ 1 and ρ 2 are two-dimensional transverse position vectors, z is the direction of propagation) and at frequency ω is given as [16]: where the matrix elements are with E α (r, ω) being a realization of the Cartesian component α of the beam-like electric field. Here * and . denote the complex conjugation and the average over the ensemble of monochromatic realizations, respectively. The single-point Stokes parameters are related to the cross-spectral density matrix elements as follows [17]: where the angular brackets are used here in contrast with a classic notation to distinguish the averaged and the instantaneous sets of parameters (the latter will be introduced below). Further, the two-point (generalized) Stokes parameters are defined as [32].
S 0 (r 1 , r 2 , ω) = W xx (r 1 , r 2 , ω) + W yy (r 1 , r 2 , ω), S 2 (r 1 , r 2 , ω) = W xy (r 1 , r 2 , ω) + W yx (r 1 , r 2 , ω), (4c) On assuming that a beam propagates in free space from the source plane z = 0 into half-space z > 0, close to the positive z-direction, and that the realizations of the electric field components in the source plane and in the receiver plane, at the point is r = (ρ, z), are E s α (ρ , ω) and E α (ρ, z, ω) respectively, the relation between them can be expressed via the Huygens-Fresnel integral [17,33]: where G(ρ − ρ , z, ω) is the free-space Green's function of a point source: with k = ω/c and c being the speed of light in a vacuum. Then the elements of the crossspectral density matrices at a pair of points in the source plane and in the receiver plane are related as: If the half-space z > 0 is filled with classic (homogeneous, isotropic, Kolmogorov-like) atmospheric turbulence, the product of the Green's functions become its second-order correlation function [6]: where Φ n (κ) is the turbulence power spectrum and κ is the spatial frequency. This implies that the cross-spectral density matrix elements become: If the source radiating the beam is of the electromagnetic Gaussian Schell model class, its cross-spectral density matrix has the form [17]: Here A α and A β are the amplitudes of the electric field vector components E α and E β , B αβ is the correlation coefficient between E α and E β components, B xy = B xy exp(−jφ), φ is the phase of the xy component, B xy = B * yx , σ 2 is the variance of the source spectral density width and δ 2 αβ is the variance of the correlation between components E α and E β . The parameters of the source must satisfy the beam conditions and the realizability conditions that are the consequence of the non-negative definiteness and quasi-Hermiticity of the cross-spectral density matrix [17]. Substituting from Equation (10) into Equation (9) yields: is the refractive-index structure constant, κ m = 5.92/l 0 , l 0 is the Kolmogorov inner scale length [34], one finds that P = π 2 k 2 z Changing variables as ρ s = 1 2 (ρ 1 + ρ 2 ) and ρ d = (ρ 1 − ρ 2 ), and after performing integrations, Equation (12) reduces to form (see Appendix A): where

The Correlations in Stokes Parameters in Turbulence
The fluctuation in each instantaneous Stokes parameter about its average value can be defined as: To cover all possible correlations among the Stokes fluctuations, the 4 × 4-dimensional correlation matrix C(r 1 , r 2 , ω) is defined [27]: with n,m = 0,1,2,3. The components of this matrix can be also expressed as [27]: where a, b, c, d = x, y, and are the Pauli spin matrices. We can find the fourth-order coherence of the field in the receiver plane (in Equation (16)) based on the extended Huygens-Fresnel integral [26]. In addition to the second-order correlations that were discussed in Section 2, the right side of Equation (16) involves the fourth-order correlation functions. The relation between these correlation functions in the source plane and in the field takes the form: where the last term gives the fourth-order correlation function of complex phase ψ induced by the refractive index fluctuations in turbulence, and subscript m denotes average taken over the ensemble of turbulence realizations.
If the source fluctuations obey the Gaussian statistics, the fourth-order correlations in the source plane can be written in terms of the second-order correlations by using the Gaussian moment theorem [5]: Under these assumptions, the correlations in Equation (19) (in the source plane) become: Further, under the conditions of weak turbulence, we can express the last exponential term in Equation (18) as [35]: where with σ 2 χ = 0.124k 7/6 C 2 n z 11/6 being the log-amplitude variance of the scintillation, ρ χS = 0.114C 2 n k 13/6 z 5/6 −1/2 and ρ 0 = 0.545C 2 n k 2 z −3/5 being the transverse coherence radius of the spherical wave. Since the on-axis scintillation is numerically examined in this study, setting ρ 1 = ρ 2 = ρ 3 = ρ 4 = 0, yields the analytic formula for the fourth-order coherence matrix (see Appendix B): where the arguments z and ω of γ n on the right side were omitted for brevity.
Thus, the on-axis Stokes scintillations can be found with Equation (16) by using the second and the fourth-order coherence functions as: For scalar fields, the intensity-intensity correlation function is conventionally normalized by the product of the average intensities of the field. However, on considering the fluctuations in the electromagnetic fields, there are several ways of introducing the normalization. For example, normalization was used in [27] where S 0 (ρ 1 , ω) and S 0 (ρ 2 , ω) are: We will follow an alternative normalization: where we use Equation (3a) to find the denominator in Equation (30), S 0 (ρ 1 , ω) S 0 (ρ 2 , ω) is found to be as: We then can find the cross-spectral density matrices in Equations (25) and (31) by using the analytical expression given in Equation (16).

Results
The following sets of parameters were used for all the numerical examples, unless other values are listed in the figure captions: Figure 1A: λ = 0.628 µm, σ = 2.5 cm, δ xx = 1.5 mm, δ yy = 2 mm, δ xy = δ yx = 2.5 mm, A x = 1.5, A y = 1, B xy = B * yx = 0.3e −iπ/6 , C 2 n = 1 × 10 −16 m −2/3 , l 0 = 1 mm.      Figures 1A and 3 except for δ xx , δ yy , δ xy and δ yx are multiplied by a coefficient. Figure 1 shows the average Stokes parameters S n (n = 0, 1, 2, 3) of the most general EM GSM beam propagating in the weak atmospheric turbulence. In Figure 1A, the value of C 2 n is fixed while the propagation distance L is set to vary. The decay rates of the Stokes parameters depend on the choice of the source r. m. s. correlations. Figure 1B shows the average Stokes parameters at a fixed distance L with varying values of the C 2 n . Starting from the value C 2 n = 10 −13 m −2/3 and higher, the average Stokes parameters decrease to zero due to the effects of turbulence. Figure 2 shows the 4 × 4 scintillation matrix elements σ 2 I,nm for a fixed range L but for varying C 2 n values. Figure 2A presents the range from C 2 n = 10 −18 m −2/3 (negligible turbulence) to C 2 n = 3 × 10 −15 m −2/3 (strong turbulence). In the interval C 2 n = 1 × 10 −15 m −2/3 to C 2 n = 3 × 10 −15 m −2/3 , the scintillation index matrix elements exhibit sharp rises. The sharp increase in the scintillations in Figure 2 occurring for the C 2 n values about C 2 n = 3 × 10 −15 m −2/3 can be explained by the joint effect of the source and the turbulence correlations when they become comparable in strength at the given distance (L = 500 m). Figure 2B gives the version of Figure 2A, zoomed to values of C 2 n at which the scintillation changes begin. While eight curves are shown to substantially deviate from zero, the other eight curves remain very close to it. Figure 3 presents the changes in the 4 × 4 matrix of the scintillation indices σ 2

I,nm
varying with propagation distance L, for the fixed value of C 2 n . In Figure 3A, the range of 0 to 10 km is included with two distinctive changes appearing in all the 16 matrix elements. The first set of changes occurs at distances smaller than 100 m. These changes are induced by source correlations and involve one, two, or three maxima for each of the 16 elements, the number depending on index nm. From 100 m to approximately 1 km, all the elements of the scintillation matrix take on very small values. However, starting in the 1 to 2 km range they start deviating from zero again, the effect attributed to the atmospheric turbulence. Figure 3B,C highlights the source-induced and the turbulenceinduced scintillations, respectively, within the corresponding ranges.  In Figure 4, the effect of the source correlation widths δ xx , δ yy , δ xy and δ yx on the changes in the elements of the 4 × 4 scintillation matrix is shown: for subfigures (A)-(E) the values of these parameters are multiplied by 5, 8, 10, 12, and 15 times as compared with those in Figure 4A, respectively. As source correlation widths grow, the effect of the source occurs at larger propagation distances from the source. If in Figure 4A the source-induced and the correlation-induced effects are still very well separated, in Figure 4B-D they become partially mixed, and in Figure 4E they mix completely.

Discussion
We will now give some general remarks on the behavior of the statistics of interest. First of all, we have limited ourselves only to a single EM GSM source parametric setup. Nevertheless, our choice of parameters has allowed us to illustrate the non-trivial behavior in all 16 entries of the scintillation matrix. The choice of other combinations of the 10 parameters of the EM GSM source and the application of different normalizations leads to quantitatively different but qualitatively similar behavior of the scintillation index matrix. In particular, in the limiting cases of unpolarized and uncorrelated sources as well as of sources with no phase lag between the two electric field components, some elements of the scintillation matrix vanish and the rest might yield monotonic behavior (not shown). We also point out that the scintillation matrix has at least one non-trivial element already at the source plane.
We also note that the behavior of the average Stokes parameters and the elements of the scintillation index matrix are qualitatively different: the former tend to zero with the growing propagation distance while the latter can saturate at non-zero values. According to Equations (3a-d) the average, one-point Stokes parameters are the linear combinations of the single-point cross-spectral density matrix (CSDM) elements. The diagonal elements are just the average intensities of the x and y electric field components and the off-diagonal elements are the single-point x-y electric field correlations. The off-diagonal elements can be put into the diagonal CSDM elements by a unitary transformation [36]. Hence, the CSDM elements and, therefore, the average Stokes parameters behave similarly to the average intensity. In free space, and more so in atmospheric turbulence, the average intensity at a single point would decrease as the beam propagates. In free space, this is due to diffraction alone and in turbulence, this is a cumulative effect of free-space diffraction and long-term turbulence-induced beam spread. On the other hand, the scintillation index matrix elements are the normalized quantities (see Equation (27)), and, hence, have the same power (fourth) in the field, in all the terms appearing in the numerator and the denominator of the expression.

Conclusions
We have formulated and derived the propagation laws for the 4 × 4 matrix consisting of the scintillation indices appearing on correlating the two electric field components of an optical field, generated by a stationary, electromagnetic source and propagated in a weak, isotropic, and homogeneous atmospheric turbulence. We have applied the formulation to the EM GSM sources and the Tatarskii power spectrum of atmospheric refractive index fluctuations. For treating propagation in turbulence, we used the extended Huygens Fresnel method.
We found that in the most general case of the EM GSM source, all the sixteen elements of the (on-axis) scintillation matrix are modulated differently. We also showed that two qualitatively different types of modulation occur: first, due to the source correlations and, second, due to the atmospheric turbulence. The turbulence-induced changes in the scintillation indices start appearing (for typical atmospheric fluctuations) at ranges on the order of kilometers and are monotonic. However, the changes caused by the source correlations are more complex. First of all, for sufficiently small source correlation widths they can be made appear at relatively small ranges and be completely decoupled from the turbulence-induced scintillations. In such cases, all the scintillation indices first attain their maxima/minima values, then return back to zero, and then undergo a second monotonic deviation from zero. If the source correlation widths are scaled up, the two scintillation changes can be made to coincide.
The analysis of this paper is instrumental for the understanding the use of the stationary electromagnetic beams for free-space optical communication systems, and possibly remote sensing and imaging systems, operating in the presence of atmospheric turbulence. The new diversity mechanism using the 16-channel optical carrier can be employed for increasing the quality and/or the data transmission rate in the presence of random environments. Similar treatment can be implemented for signal transmission through underwater turbulent environments or in soft, homogeneous and isotropic biological tissues.