1. INTRODUCTION

This paper lays out a methodology for solving problems of electromagnetic wave propagation in turbulence. Until recently, we (the authors) would have had two choices to find a quantity like jitter variance: apply an existing analytic result from the literature or model the propagation numerically and generate statistics through Monte Carlo simulations. While a number of such equations exist in the literature [1–3], and Sasiela formalized existing approaches into a general-purpose methodology for generating such expressions [4], we have found that a great deal of mathematical sophistication is necessary to apply those methods to anything novel. We have formulated a different approach to deriving these sorts of analytic expressions, built on equations and concepts that will be familiar to any reader with a background in paraxial wave-optics modeling. We lay out that approach in this paper, providing a method to understand and generate analytic expressions for the statistics of atmospheric propagation.

The central concept of our approach is that, just as we build up atmospheric turbulence in wave optics by summing appropriately weighted Fourier components, we can calculate atmospheric effects by first understanding the effect of a single sinusoid and then summing the appropriately weighted results. The development in this paper has many commonalities with that found in Sasiela, and it was in fact a study of that text that led to our current understanding. The advantage of this approach over previous treatments is that it will allow a wider audience that is not comfortable with the methods of random field theory and stochastic calculus to solve problems of practical interest; readers who are comfortable with existing analytic methods may also find that additional insight can be gained from our approach. We will draw on equations, literature, and phenomena likely to be familiar to those who study atmospheric propagation, such as the paraxial propagator and the methods used to generate turbulence phase screens, with each step of the development anchored in physical intuition. The methods described in this paper do involve approximations that become inaccurate for scenarios with Rytov (log-amplitude) variance $\gtrsim 0.3$, which is the same limit faced by analytic techniques derived from the first-order Rytov approximation [5].

The primary purpose of this paper, rather than to derive new results, is to explain our new methodology and to demonstrate its validity via agreement with well-established results. We will work through three examples of increasing complexity, laying the groundwork for important quantities and finding several expressions that can be compared with results in the literature. We will begin with gradient tilt, using light from a point-source beacon as an example case while developing the general methods and master equations for finding variances and power spectral densities (PSDs). We then work with Zernike tilt, using two offset beams as an example while developing the general methods to handle anisoplanatic effects. Finally, we will show how to derive expressions for the centroid motion of a generic beam using centroid tilt, arriving at results that do not currently exist in the literature. Using the equations and methods developed in our examples, readers who were daunted by previous methods will have new tools to derive analytic variances and PSDs for novel problems applicable to their research.

2. GRADIENT TILT VARIANCE AND POWER SPECTRUM

In this section, we take as an example light from a point-source beacon that travels a distance $L$ through atmospheric turbulence and passes through a circular aperture. We then wish to find the variance and power spectral density of the gradient tilt (G tilt) of the field in that aperture. In wave optics, we would calculate this by setting up the initial point-source field, propagating it through a series of discrete phase screens, clipping and collimating at the aperture, and then calculating G tilt on the field; this path is shown schematically in Fig. 1. The resulting time series of G tilt would give us a variance and a PSD, and by running many independent realizations of the atmosphere we would build up ensemble-averaged statistics. In this section, we will work through the mathematics to find analytic expressions for those same ensemble-averaged variances and PSDs.

 

Fig. 1. In this figure, light from a point-source beacon traces out a cone as it propagates a distance $L$ through the atmosphere to an aperture. The atmospheric turbulence is indicated schematically as a series of phase screens as we would model it in wave-optics code.

Download Full Size | PDF

A. Propagation of a Sinusoidal Phase Profile

In wave optics, we are accustomed to thinking of the effects of atmospheric turbulence as being approximated by a series of discrete phase screens as we have drawn in Fig. 1. Let us take this concept one step further by considering the effect of a single sinusoidal phase profile $\varphi$ with transverse wave vector $\mathit{\kappa }$, applied to our field at some position $z$ along the path as follows:

B. G Tilt from a Sinusoidal Phase Profile

G tilt is proportional to the average wavefront gradient measured over some region $\mathcal{R}$. For a generic complex field $\psi (\mathit{r})=a(\mathit{r})\exp [i\varphi (\mathit{r})]$, the component of G tilt in the $\widehat{\mathit{m}}$direction is given by

C. G-Tilt Variance from a Turbulent Phase Profile

Now that we have an expression for the G-tilt variance measured at distance $L$ due to a single sinusoidal phase profile at $z$, the next step is to use these sinusoids to construct a turbulent phase profile. From wave-optics modeling, this process is quite familiar: we model turbulence by adding together Fourier components with random phase offsets and amplitudes scaled by the appropriate power spectrum. Since G tilt is found by a linear operator acting on the phase, the G tilt due to this weighted sum of sinusoids is a sum of the G tilts caused by those sinusoids, as given in Eq. (5), with the same weights. Similarly, the fact that these sinusoids have independent random phases means that the ensemble-averaged G-tilt variance is the sum of the variances contributed by each sinusoid, as calculated with Eq. (6), again with the same spectral weighting as used to build the phase screen from Fourier components.

In Appendix B we show in detail how the methods that are used in wave optics for generating phase screens can be translated into the integral that is needed for this analytic approach. The key result for G-tilt variance, found in Eq. (B7), is

D. Path-Integrated G-Tilt Variance

The final step in our derivation is to integrate contributions to the G-tilt variance at $z=L$ from all turbulent phase profiles along the path of propagation. Mathematically, this step is trivial: integrating Eq. (7) from $z=0$ to $z=L$ and pulling terms as far out of the nested integrals as possible, we obtain

E. G-Tilt Power Spectral Density

Once we have an expression for the variance of a quantity of interest like G tilt, we can obtain an expression for its PSD with very little additional effort: the conversion of a variance expression to a PSD expression can be accomplished with a change of variables. Other related techniques exist, such as finding a covariance expression and Fourier transforming it [1], but the method described here allows us to arrive at the required mathematical transformation directly from our intuition about wave-optics modeling.

The first thing that we need to do is introduce time dependence into our phase profiles. We adopt Taylor’s frozen turbulence hypothesis, which allows us to assume that all temporal variation is attributable to turbulent phase drifting through the path of propagation at apparent wind velocity $\mathit{v}(z)$. Mathematically, as we would do in wave-optics modeling, this means that we add a time-dependent phase shift to the sinusoid in Eq. (1), $\alpha (\mathit{\kappa },z)\to \alpha (\mathit{\kappa },z)-\mathit{\kappa }·\mathit{v}(z)t$. This substitution for $\alpha$ carries through our propagation calculations, so that the phase profile after propagation from $z$ to $L$ is

3. MASTER EQUATIONS FOR VARIANCES AND PSDS

In the previous section, we showed how to derive expressions for the variance and PSD of G tilt integrated over an atmospheric channel. Though derived for G tilt, the expressions we found are trivially generalizable into master equations for variances and PSDs of other quantities. Dropping the subscripts from Eq. (9) gives us the variance master equation

Similarly, Eq. (20) becomes the PSD master equation

To apply these master equations to any new quantity, we simply need to find the variance ${\sigma }^2(\mathit{\kappa },z;L)$ that we measure at $L$ due to a single sinusoidal phase profile at $z$. We have already seen how to do this for G tilt; the following sections will show how to do it for Zernike tilt (Z tilt) or other Zernike modes, anisoplanatic differences between quantities, C tilt, and centroid jitter.

4. ZERNIKE TILT ANISOPLANATISM

In Section 2, we showed how to derive expressions for the variance and PSD of a quantity integrated over an atmospheric channel using G tilt as an example. A similar derivation for Z tilt would use the master equations in Section 3 and follow the G-tilt development so closely that there is not much additional value in going through it here. There are, however, many quantities that we might be interested in calculating that are combinations of simpler measurements; in this section, we will increase complexity, showing the method for finding anisoplanatic effects. As our example, we will find the difference in Z tilt for two beams propagating along paths with a spatial offset as illustrated in Fig. 2.

 

Fig. 2. For Z-tilt anisoplanatism, we consider two beams that propagate along paths displaced from one another in space. In this illustration, we have two point sources propagating to a common aperture. Both beams follow cones (solid lines) with the same $z$-dependent diameter $D(z)$, with center lines (dashed lines) separated by the distance $d(z)$. Because the two cones do not completely overlap, the Z tilts picked up by the two beams on the way to the aperture will differ, with the difference being the Z-tilt anisoplanatism.

Download Full Size | PDF

A. Z Tilt from a Sinusoidal Phase Profile

Since we have not worked with Z tilt yet in this paper, we will begin by finding the basic expressions for the Z tilt and Z-tilt variance of a single beam. Z tilt is the amount of overall linear phase tilt present in a field over a circular aperture; as with G tilt, this is a linear, phase-only measurement. Our first step, as with G tilt, is to determine the Z-tilt contribution from a sinusoidal phase profile and then determine the variance in the tilt from the ensemble of potential phases for that profile. We will again find the tilt in the direction of a unit vector $\widehat{\mathit{m}}$ oriented with the angle ${\phi }_m$.

We measure the Z tilt on the beam by making use of the orthogonality of the Zernike modes on the circle, so we can determine the Z tilt present in a phase profile by taking the overlap of the phase with a simple phase ramp

As we did with G tilt, we can now turn this into an ensemble-averaged variance by calculating the average value of ${\theta }_{\mathrm{Z},m}^2$ over all possible values of the phase $\alpha$ as follows:

B. Differential Z Tilt from a Sinusoidal Phase Profile

To get anisoplanatic effects, we compare beams with offsets either in space or time. In wave-optics modeling, this would mean propagating the two beams through the same set of phase screens, where they would pick up different Z tilts due to the offset. Analytically, any shift of the field in space or in time shows up in the equations as a modification to the phase $\alpha (\mathit{\kappa },z)$, again because the beam passes through a different portion of the sinusoidal phase profile. For a spatial offset, we take $\alpha (\mathit{\kappa },z)\to \alpha (\mathit{\kappa },z)+\mathit{\kappa }·\mathit{d}(z)$ for the initial phase in Eq. (1). This altered $\alpha$ term is unchanged by propagation to $L$ and the calculation of Z tilt, so we make the same substitution in Eq. (27) and have

C. Path-Integrated Z-Tilt Anisoplanatism

In Eq. (33) we have the expression for the differential Z-tilt variance from a sinusoidal profile for two beams with an offset $\mathit{d}(z)$. With this in hand, calculating a path-integrated variance and PSD is a straightforward application of the master equations in Section 3; the fact that this single-sinusoid variance is for an anisoplanatic effect does not matter.

By making appropriate choices for $D(z)$ and $d(z)$ for the beams, we can derive a number of different anisoplanatic expressions. For example, as shown in Fig. 2, taking $D(z)=(z/L)D$ and $d(z)=(1-z/L)d$ allows us to find the differential jitter between two point sources separated by a distance $d$ but measured in the same aperture. At this stage, however, we do not need to specify these forms before we derive the path-integrated expression through turbulence. We take Eq. (22) and substitute Eq. (33) for ${\sigma }^2(\mathit{\kappa },z;L)$ and Eq. (8) for ${\mathrm{\Phi }}_n^0(\kappa )$ to have

5. CENTROID TILT AND SCORING BEAM JITTER VARIANCE

Centroid motion on a target is directly related to centroid tilt (C tilt) and not the G or Z tilts typically considered in the literature. Because jitter of a directed-energy beam in the target plane is a central quantity in the analysis of directed-energy applications, we feel it important to develop the corrected expressions in this section. We do not know whether previous jitter analyses relied on Z tilt ([4], Section 4.7) because of an incorrect assumption that Z tilt is the quantity directly related to centroid motion or because the previous analytical approaches allowed for the treatment of Z tilt but not C tilt.

Unlike G and Z tilts, both of which are phase-only measurements, C tilt includes amplitude information, meaning that some of our initial treatment of it will differ from the previous sections. C tilt can be defined multiple ways; the most intuitive definition is that the centroid of the beam moves in the direction given by the C tilt, and C tilt gives the location of the centroid if that beam is propagated to focus. The most mathematically relevant definition for our development is that it is the irradiance-weighted gradient of the beam in the $\widehat{\mathit{m}}$ direction

A. C Tilt from a Sinusoidal Phase Profile

For G and Z tilts, we started by finding the phase profile at $L$, where we were measuring the tilt, that resulted from a sinusoidal profile at $z$. That phase profile and its tilt content were changed by propagation from $z$ to $L$, and tilt had to be measured with the phase as it appeared in the aperture. Since G and Z tilts are phase-only quantities, this analysis sufficed. C tilt, however, includes both phase and amplitude information, so we would need to ask not only how the sinusoidal phase is changed by propagation but also how that phase converts to amplitude and how $a^2(\mathit{r})$ itself changes. We are saved from this more complex analysis, however, by the fact that C tilt is unaffected by propagation in the paraxial limit (proven in Appendix C.B). This means that we can do our C-tilt calculation at $z$ without having to know how the field changes in propagating to $L$.

Applying Eq. (37) to the sinusoidal phase profile from Eq. (1), we have

B. Centroid Motion

In the previous section, we found expressions for the C tilt imparted to a beam by a sinusoidal phase profile. C tilt indicates the direction that the centroid is moving; the actual measured motion of the centroid after propagation is this tilt times the propagation distance. Thus, the centroid motion measured at $L$ in direction $\widehat{\mathit{m}}$ due to the sinusoidal phase at a position $z$ is

C. Centroid Jitter Variance on a Target

We are now in a position to find the quantity that motivated interest in C tilt to begin with, the centroid jitter variance of a beam in the target plane. We take Eq. (22) for our path-integrated variance and substitute Eq. (43) for ${\sigma }^2(\mathit{\kappa },z;L)$ and the modified von Karman ${\mathrm{\Phi }}_n^0(\kappa )$ in Eq. (8) to have

D. Gaussian-Beam Centroid Jitter

For the special case of a circular Gaussian beam, the irradiance profile has the simple analytic form

If we take Eq. (47) in the limit of Kolmogorov turbulence, the inner integral has an analytic solution, leaving us with

To test these new analytic results, we have made comparisons between them and Monte Carlo wave-optics results for beam jitter, shown in Fig. 3. For a top-hat beam, we get good agreement with our analytic results in Eq. (45) up through a spherical-wave Rytov variance of 0.3. Between a Rytov variance of 0.3 and 0.4, the analytic curves begin to differ from the wave-optics results as the weak-turbulence approximations in our analytic treatment break down. Figure 3 also shows results for a focused Gaussian beam, using the known analytic expression for $w(z)$ and taking $w(0)=2^{-3/2}D$ to approximate the far field of the top-hat beam. We again see good agreement between wave optics and the analytic results here, using Eq. (48), with disagreement appearing at the same strong-turbulence regime. The top-hat analytic results are a constant factor of 0.874 times the Gaussian analytic results over the full range of the plot, suggesting that the much simpler result in Eq. (48) can be used as an approximation for the centroid jitter of other beam shapes with modest rescaling.

 

Fig. 3. Results for centroid jitter variance on target for focused beams that have propagated through different strengths of atmospheric turbulence. In all cases we took a 1 μm wavelength beam and propagated 3 km through Kolmogorov turbulence with constant $C_n^2$. For the top-hat beam we used a 20-cm aperture, and for the Gaussian we used $w=2^{-3/2}\times 20\mathrm{cm}$ to match the top hat. The red and purple theory curves lie largely on top of their corresponding wave-optics results until we violate the weak-turbulence approximations in the theory for Rytov variance $\gtrsim 0.3$.

Download Full Size | PDF

6. CONCLUSIONS

In this paper, we have worked through the derivations for calculating G, Z, and C tilts, the centroid motion resulting from C tilt, and anisoplanatic effects. In all cases, we started with the effects of simple sinusoidal phase profiles and built up to path-integrated variances or PSDs. We have attempted to base that development on physical intuition or common wave-optics techniques to make it more accessible to readers with backgrounds like our own. The methods here and master equations in Section 3 can be applied for any phenomenon of interest, as long as that quantity is linear in the phase and can be calculated for a generic sinusoidal phase profile.

Though this approach is based on different starting points and intuition, the terms we have found are analogous to those in previous work, and our results are identical when calculated for the same conditions. We have, however, been able to derive expressions for C tilt on a generic beam, which were not present in previous analytic treatments and which allow us to find the correct centroid motion of a beam. There are many more results that could be found with the methods described here, and it is our hope that this approach will allow others to pursue such work in the future, starting with different intuition and backgrounds than assumed by other approaches in the literature.

APPENDIX A: THE PROPAGATION OF SINUSOIDAL PHASE PROFILES

In order to calculate some effects due to sinusoidal phase profiles, we need to know how those profiles change due to propagation from where they are applied to where they are measured. There are three classes of fields of interest: diverging spherical waves, converging spherical waves, and plane waves. We begin at $z=z_0$ with a field

(A1)$$\psi (\mathit{r},z_0)=\exp \left(i\frac{k}{2f_0}{r}^2\right)\exp [iϵ\sin (\mathit{\kappa }·\mathit{r}+\alpha )],$$

where the parabolic phase of focal length $f_0$ in the first term can cover all three cases of interest, and the second term has our sinusoidal phase profile with amplitude $ϵ$, wave vector $\mathit{\kappa }$, and phase offset $\alpha$. We have defined the sign of $f_0$ such that a diverging beam, as illustrated in Fig. 4, has $f_0>0$.

Using the Jacobi–Anger expansion, we express the field in Eq. (A1) as the infinite sum

(A2)$$\psi (\mathit{r},z_0)=\exp \left(i\frac{k}{2f_0}{r}^2\right)\sum _{n=-\infty }^{\infty }{J}_n(ϵ)e^{in\mathit{\kappa }·\mathit{r}}{e}^{in\alpha }.$$

Because propagation is a linear operation, we can separately consider the propagation of each component in this sum. Applying the paraxial propagation equation {[7], Eq. (4-14)}, we have

(A3)$$\begin{gathered} \psi (\mathit{r},z_0+z)=\frac{e^{ikz}}{i\lambda z}\sum _{n=-\infty }^{\infty }{J}_n(ϵ)e^{in\alpha }\iint \mathrm{d}{\mathit{r}}_0{e}^{in\mathit{\kappa }·{\mathit{r}}_0}\exp \left(i\frac{k}{2f_0}{r}_0^2\right) \\ \times \exp \left[i\frac{k}{2z}{({\mathit{r}}_0-\mathit{r})}^2\right]. \end{gathered}$$

After propagating a distance $z$, we could expect the field to have a new focal length $f_1=f_0+z$ and to have changed in size by a transverse factor $s=f_1/f_0=D(z_0+z)/D(z_0)$ based on the simple ray-optics intuition suggested by Fig. 4. Using these definitions of $f_1$ and $s$, it can be shown algebraically that

(A4)$$\frac{k}{2f_0}{r}_0^2+\frac{k}{2z}{({\mathit{r}}_0-\mathit{r})}^2=\frac{k}{2f_1}{r}^2+\frac{ks}{2z}{({\mathit{r}}_0-\frac{\mathit{r}}{s})}^2.$$

Using this substitution in Eq. (A3), we then have

(A5)$$\begin{gathered} \psi (\mathit{r},z_0+z)=\frac{e^{ikz}}{i\lambda z}\exp \left(i\frac{k}{2f_1}{r}^2\right)\sum _{n=-\infty }^{\infty }{J}_n(ϵ)e^{in\alpha } \\ \times \iint \mathrm{d}{\mathit{r}}_0{e}^{in\mathit{\kappa }·{\mathit{r}}_0}\exp \left[i\frac{ks}{2z}{({\mathit{r}}_0-\frac{\mathit{r}}{s})}^2\right]. \end{gathered}$$

In this form, the double integral is now the Fourier transform of the paraxial propagator with some variable changes and rescaling. Performing that Fourier transform ([7], Table 2.1) gives us

(A6)$$\begin{gathered} \psi (\mathit{r},z_0+z)=\frac{e^{ikz}}{s}\exp \left(i\frac{k}{2f_1}{r}^2\right)\sum _{n=-\infty }^{\infty }{J}_n(ϵ)e^{in\mathit{\kappa }·\mathit{r}/s}{e}^{in\alpha } \\ \times \exp (-i\frac{n^2{\kappa }^{2}z}{2ks}). \end{gathered}$$

This is the exact expression for our propagated field. To make further analytic progress, we need to make a small-$ϵ$ approximation to this result. Most operations in this paper are strictly linear in the phase, relying at most on the approximation of paraxial propagation. Our conclusions here about how sinusoidal phase profiles propagate are the exception where we need to make a weak-turbulence approximation to ensure linearity. Because our G- and Z-tilt calculations rely on these results, that makes all of them dependent on that approximation. In this limit we have $J_n(ϵ)\propto {ϵ}^{|n|}$, and the series to the first order in $ϵ$ is equal to the first-order expansion of the more useful form

(A7)$$\begin{gathered} \psi (\mathit{r},z_0+z)\approx \frac{e^{ikz}}{s}\exp \left(i\frac{k}{2f_1}{r}^2\right) \\ \times \exp \left[iϵ\cos \left(\frac{{\kappa }^{2}z}{2ks}\right)\sin (\frac{\mathit{\kappa }·\mathit{r}}{s}+\alpha )\right] \\ \times \exp \left[ϵ\sin \right(\frac{{\kappa }^{2}z}{2ks}\left)\sin (\frac{\mathit{\kappa }·\mathit{r}}{s}+\alpha )\right]. \end{gathered}$$

The first and second complex exponentials in this result are the global piston and focus after propagation, respectively. The third complex exponential gives the new sinusoidal phase profile. If we identify the normalized phase $\varphi (\mathit{\kappa },\mathit{r})=\sin (\mathit{\kappa }·\mathit{r}+\alpha )$ in Eq. (A1), we then have

(A8)$$\begin{gathered} \varphi (\mathit{\kappa },\mathit{r},z_0;z_0+z)=\sin (\frac{D(z_0)}{D(z_0+z)}\mathit{\kappa }·\mathit{r}+\alpha ) \\ \times \cos \left(\frac{{\kappa }^{2}z}{2k}\frac{D(z_0)}{D(z_0+z)}\right) \end{gathered}$$

in Eq. (A7), where we have denoted that the phase is applied at $z_0$ but measured after propagation to $z_0+z$, and we have substituted in one of the definitions of $s$. We see that our original sinusoidal phase profile has been transversely rescaled and reduced in amplitude by the cosine factor. The final exponential in Eq. (A7) shows that this loss in amplitude of the phase profile is matched by a gain in a log-amplitude profile, such that the sum of phase variance and log-amplitude variance is conserved. The transverse rescaling of our phase profile is easily understood by the ray-optics picture in Fig. 4. The change in amplitude of the phase profile can be related to a known diffraction phenomenon if we rewrite the cosine amplitude factor in Eq. (A8) as

(A9)$$\cos \left(\frac{{\kappa }^{2}z}{2k}\frac{D(z_0)}{D(z_0+z)}\right)=\cos \left(2\pi \frac{z}{z_{\mathrm{T}}}\frac{D(z_0)}{D(z_0+z)}\right),$$

where $z_{\mathrm{T}}=4\pi k/{\kappa }^2$ is the Talbot length. The well-known Talbot effect [7–10] occurs when a field with periodic structure sees that structure change under propagation but return to an exact copy of the original field after a propagation distance $z_{\mathrm{T}}$ ([7], Section 4.5.3). We have shown that it holds for spherical wave fronts, provided that the Talbot length is modified by a geometric rescaling factor.

 

Fig. 4. Illustration of our setup to calculate the effect of propagation on a sinusoidal phase profile. We start with a field of focal length $f_0$, shown here as a diverging beam, with a sinusoidal phase profile on top of the parabolic phase. In propagating a distance $z$, our initial phase $\varphi (\mathit{\kappa },\mathit{r})$ changes to $\varphi (\mathit{\kappa },\mathit{r},z_0;z_0+z)$. We will derive an expression for the changes to this phase profile from propagation.

Download Full Size | PDF

APPENDIX B: BUILDING ATMOSPHERIC TURBULENCE FROM SINUSOIDS

The development of propagation statistics in this paper starts from the calculation of measurable effects from simple sinusoidal phase profiles. To apply these results to atmospheric turbulence, we need expressions that allow the combination of sinusoidal profiles as Fourier components of atmospheric turbulence with appropriate weighting for the turbulence spectrum.

We begin with the equation for generating a random realization of the phase change $\mathrm{\Delta }\varphi$ picked up by a propagating field due to a slab of atmosphere of thickness $\mathrm{\Delta }z$ at location $z$, given by [11]

(B1)$$\mathrm{\Delta }\varphi (\mathit{r},z)=\iint \mathrm{d}\mathit{\kappa }\frac{g(\mathit{\kappa },z)}{q^{1/2}}{[F_{\mathrm{\Delta }\varphi }(\kappa ,z)]}^{1/2}{e}^{i\mathit{\kappa }·\mathit{r}},$$

where $g(\mathit{\kappa },z)$ is a realization of delta-correlated, complex, Hermitian, Gaussian white noise with mean-square value $q$, and $F_{\mathrm{\Delta }\varphi }$ is the power spectral density of $\mathrm{\Delta }\varphi$, which we have assumed to be isotropic. If we let $B(\mathit{\kappa },z)=|g(\mathit{\kappa },z)/q^{1/2}|$ and $\alpha (\mathit{\kappa },z)=\arg [g(\mathit{\kappa },z)]+\pi /2$, and we integrate over the half-plane to handle the Hermitian symmetry, we have

(B2)$$\begin{gathered} \mathrm{\Delta }\varphi (\mathit{r},z)=2{\int }_{-\infty }^{\infty }\mathrm{d}{\kappa }_x{\int }_0^{\infty }\mathrm{d}{\kappa }_{y}B(\mathit{\kappa },z){[F_{\mathrm{\Delta }\varphi }(\kappa ,z)]}^{1/2} \\ \times \sin [\mathit{\kappa }·\mathit{r}+\alpha (\mathit{\kappa },z)]. \end{gathered}$$

Our next steps depend only on measuring a quantity that is linear in the phase; taking G tilt as an example, we find that the G tilt applied by this realization of atmospheric turbulence, if applied at $z$ and measured at $L$, is given by

(B3)$$\begin{gathered} {\theta }_{\mathrm{G},m}(z;L)=2{\int }_{-\infty }^{\infty }\mathrm{d}{\kappa }_x{\int }_0^{\infty }\mathrm{d}{\kappa }_{y}B(\mathit{\kappa },z){[F_{\mathrm{\Delta }\varphi }(\kappa ,z)]}^{1/2} \\ \times \cos [\alpha (\mathit{\kappa },z)]{\tilde{\theta }}_{\mathrm{G},m}(\mathit{\kappa },z;L), \end{gathered}$$

where we have introduced ${\theta }_{\mathrm{G},m}=\cos [\alpha (\mathit{\kappa },z)]{\tilde{\theta }}_{\mathrm{G},m}(\mathit{\kappa },z;L)$ so that ${\tilde{\theta }}_{\mathrm{G},m}$ contains all of the terms in ${\theta }_{\mathrm{G},m}$ that do not vary with $\alpha$ [see Eq. (5)].

We now need to obtain the ensemble-averaged G-tilt variance ${\sigma }_{\mathrm{G},m}^2(z;L)=⟨{\theta }_{\mathrm{G},m}^2(z;L)⟩\mathrm{.}$ Because $B(\mathit{\kappa },z)$ and $\cos [\alpha (\mathit{\kappa },z)]$ are independent at every $\mathit{\kappa }$ in the half-plane, the ensemble averages can be brought inside the integrals and evaluated to give us

(B4)$${\sigma }_{\mathrm{G},m}^2(z;L)=4{\int }_{-\infty }^{\infty }\mathrm{d}{\kappa }_x{\int }_0^{\infty }\mathrm{d}{\kappa }_y{F}_{\mathrm{\Delta }\varphi }(\kappa ,z){\sigma }_{\mathrm{G},m}^2(\mathit{\kappa },z;L),$$

where ${\sigma }_{\mathrm{G},m}^2(\mathit{\kappa },z;L)$ is the G-tilt variance for a single sinusoidal profile over an ensemble of phase values as found in Eq. (6). Restoring the integral to be over the full plane gives us

(B5)$${\sigma }_{\mathrm{G},m}^2(z;L)=2\iint \mathrm{d}\mathit{\kappa }{F}_{\mathrm{\Delta }\varphi }(\kappa ,z){\sigma }_{\mathrm{G},m}^2(\mathit{\kappa },z;L).$$

We now introduce the form of $F_{\mathrm{\Delta }\varphi }$ for atmospheric turbulence [11],

(B6)$$F_{\mathrm{\Delta }\varphi }(\kappa ,z)=2\pi k^2{\mathrm{\Phi }}_n^0(\kappa ){\int }_z^{z+\mathrm{\Delta }z}\mathrm{d}{z}^{\prime }{C}_n^2(z^{\prime })$$

for a slab of atmosphere $\mathrm{\Delta }z$ thick at some point $z$ along the propagation path, where ${\mathrm{\Phi }}_n^0$ is the scaled refractive index power spectral density [given for the modified von Karman spectrum in Eq. (8)]. If we take the limit $\mathrm{\Delta }z\to \mathrm{d}z$ and substitute Eq. (B6) into Eq. (B5), we get the differential G-tilt contribution for an infinitesimal portion of the path at $z$ as follows:

(B7)$$\mathrm{d}{\sigma }_{\mathrm{G},m}^2(z;L)=4\pi k^2\mathrm{d}z{C}_n^2(z)\iint \mathrm{d}\mathit{\kappa }{\mathrm{\Phi }}_n^0(\kappa ){\sigma }_{\mathrm{G},m}^2(\mathit{\kappa },z;L).$$

In order to integrate both sides of Eq. (B7), we will need to neglect any atmospheric turbulence correlations along the $z$ axis. This assumption is justified as long as the propagation path spans many atmospheric correlation lengths ([12], Section 8.6). With this limitation in mind, we obtain

(B8)$${\sigma }_{\mathrm{G},m}^2=4\pi k^2{\int }_0^L\mathrm{d}z{C}_n^2(z){\int }_0^{\infty }\mathrm{d}\kappa \kappa {\mathrm{\Phi }}_n^0(\kappa ){\int }_0^{2\pi }\mathrm{d}{\phi }_{\kappa }{\sigma }_{\mathrm{G},m}^2(\mathit{\kappa },z;L),$$

which is the ensemble-averaged variance of the G tilt for propagation through atmospheric turbulence.

APPENDIX C: PROPERTIES OF CENTROID TILT

In the existing literature on atmospheric propagation statistics, gradient and Zernike tilts (G and Z tilts) are the common quantities studied. Because centroid tilt (C tilt) is less familiar, it is worth working through the derivations of a number of fundamental properties that we rely on in this work.

A. Equivalent C-Tilt Definitions

There are two definitions that we will show to be equivalent in this section: C tilt can be defined as the irradiance-weighted average phase gradient or as the tilt angle measured by bringing a beam to focus and measuring the irradiance centroid location.

We start with an arbitrary complex field $\psi (\mathit{r})$, with a total power given by

(C1)$$P=\iint \mathrm{d}\mathit{r}{|\psi (\mathit{r})|}^2=\frac{1}{4{\pi }^2}\iint \mathrm{d}\mathit{\kappa }{|\mathcal{F}[\psi (\mathit{r})]|}^2,$$

where $\mathcal{F}$ denotes a Fourier transform. Keeping in mind that, in the paraxial regime, taking the Fourier transform is equivalent to propagating to focus, we can get our centroid measurement at focus by taking

(C2)$${\mathit{X}}_{\kappa }=\frac{1}{4{\pi }^{2}P}\iint \mathrm{d}\mathit{\kappa }\mathit{\kappa }{|\mathcal{F}[\psi (\mathit{r})]|}^2.$$

Applying the identity $\mathcal{F}[\nabla \psi (\mathit{r})]=i\mathit{\kappa }\mathcal{F}[\psi (\mathit{r})]$ and dividing by $k$ to convert to an angular quantity for C tilt, we have

(C3)$${\mathit{\theta }}_{\mathrm{C}}=\frac{{\mathbf{X}}_{\kappa }}{k}=\frac{1}{i4{\pi }^{2}Pk}\iint \mathrm{d}\mathit{\kappa }\mathcal{F}[\nabla \psi (\mathit{r})]{\{\mathcal{F}[\psi (\mathit{r})]\}}^{*}.$$

By writing out the Fourier transform integrals explicitly and carrying out the leading $\mathit{\kappa }$ integral, we generate a delta function and reduce this to

(C4)$${\mathit{\theta }}_{\mathrm{C}}=\frac{1}{iPk}\iint \mathrm{d}\mathit{r}[\nabla \psi (\mathit{r})]{\psi }^{*}(\mathit{r}).$$

Substituting $a(\mathit{r})e^{i\varphi (\mathit{r})}$ for $\psi (\mathit{r})$ gives us

(C5)$${\mathit{\theta }}_{\mathrm{C}}=\frac{1}{2iPk}\iint \mathrm{d}\mathit{r}\nabla [a^2(\mathit{r})]+\frac{1}{Pk}\iint \mathrm{d}\mathit{r}{a}^2(\mathit{r})\nabla \varphi (\mathit{r}).$$

For a differentiable $a^2$ with finite power, which should be true of any irradiance profile, Green’s theorem can be used to show that the first integral above goes to zero. This leaves

(C6)$${\mathit{\theta }}_{\mathrm{C}}=\frac{1}{Pk}\iint \mathrm{d}\mathit{r}{a}^2(\mathit{r})\nabla \varphi (\mathit{r}),$$

giving the C tilt as the irradiance-weighted average phase gradient. Because C tilt is linear in the phase, as shown by this result, adding phase to a field merely adds the C tilt associated with the phase to the C tilt already on the field.

B. Effect of Propagation on C Tilt

Suppose that we have a field $\psi (\mathit{r},0)$ and propagate it in vacuum a distance $z$ to $\psi (\mathit{r},z)$. In this section we will show that the C tilts measured in these two planes are identical, without saying anything about how the irradiance profile may have changed. This will show that the paraxial propagator itself does not change C tilt on a beam, unlike G and Z tilts, both of which change under propagation.

Let us begin by considering the propagation between the planes as the convolution of the original $\psi (\mathit{r},0)$ with a propagator ([7], Section 4.2). This is most easily analyzed in Fourier space, where the convolution becomes a product

(C7)$$\mathcal{F}[\psi (\mathit{r},z)]=\mathcal{F}[\psi (\mathit{r},0)]H(\mathit{\kappa },z).$$

The transfer function $H$, accepting the validity of the Fresnel approximation, is given by {[7], Eq. (4-20)}

(C8)$$H(\mathit{\kappa },z)=\exp [i\frac{2\pi }{\lambda }z-i\frac{\lambda z}{4\pi }({\kappa }_x^2+{\kappa }_y^2)].$$

This means that

(C9)$${|\mathcal{F}[\psi (\mathit{r},z)]|}^2={|\mathcal{F}[\psi (\mathit{r},0)]|}^2.$$

Thus, the Fourier transforms of the fields before and after propagation have identical absolute squares, regardless of how the field itself may have changed. Since one definition of C tilt is the centroid of the absolute square of the Fourier transform, the fields before and after vacuum propagation have identical C tilts.

C. Effect of Propagation on Centroid Measurements

In this section we show that the propagation of an arbitrary field with a given C tilt causes the irradiance centroid to shift by the C tilt times the propagation distance. This is for arbitrary vacuum propagation, not just propagation to focus, and is a key step to allowing integration of C tilts along a path to get the centroid motion on a target.

We start with the paraxial propagator from Goodman {[7], Eq. (4-17)},

(C10)$$\psi (\mathit{r},z)=\frac{e^{ikz}}{i\lambda z}\exp \left(i\frac{k}{2z}{r}^2\right)\mathcal{F}[\psi ({\mathit{r}}_0,0)\exp \left(i\frac{k}{2z}{r}_0^2\right)],$$

where we can equate $\mathit{\kappa }=2\pi \mathit{r}/\lambda z$ after the transform. Our irradiance centroid after propagation is now given by

(C11)$$\mathit{X}(z)=\frac{1}{P}\iint \mathrm{d}\mathit{r}\mathit{r}{|\frac{e^{ikz}}{i\lambda z}\exp \left(i\frac{k}{2z}{r}^2\right)\mathcal{F}[\psi ({\mathit{r}}_0,0)\exp \left(i\frac{k}{2z}{r}_0^2\right)\left]\right|}^2.$$

Using our variable identity to transform from $\mathit{r}$ to $\mathit{\kappa }$ in the integrals, we have

(C12)$$\begin{gathered} \mathit{X}(z)=\frac{\lambda z}{8{\pi }^{3}P}\iint \mathrm{d}\mathit{\kappa }\mathit{\kappa }\mathcal{F}[\psi ({\mathit{r}}_0,0)\exp \left(i\frac{k}{2z}{r}_0^2\right)] \\ \times {\mathcal{F}}^{*}[\psi ({\mathit{r}}_1,0)\exp \left(i\frac{k}{2z}{r}_1^2\right)]. \end{gathered}$$

With our Fourier transform identity $\mathcal{F}[\nabla \psi (\mathit{r})]=i\mathit{\kappa }\mathcal{F}[\psi (\mathit{r})]$, this becomes

(C13)$$\begin{gathered} \mathit{X}(z)=\frac{\lambda z}{i8{\pi }^{3}P}\iint \mathrm{d}\mathit{\kappa }\mathcal{F}\{\nabla [\psi ({\mathit{r}}_0,0)\exp \left(i\frac{k}{2z}{r}_0^2\right)\left]\right\} \\ \times {\mathcal{F}}^{*}[\psi ({\mathbf{r}}_1,0)\exp \left(i\frac{k}{2z}{r}_1^2\right)]. \end{gathered}$$

If we write out the Fourier integrals explicitly and integrate over $\mathit{\kappa }$, we get a delta function and can reduce this to

(C14)$$\begin{gathered} \mathit{X}(z)=\frac{\lambda z}{i2\pi P}\iint \mathrm{d}\mathit{r}\nabla [\psi (\mathit{r},0)\exp \left(i\frac{k}{2z}{r}^2\right)] \\ \times {\psi }^{*}(\mathit{r},0)\exp (-i\frac{k}{2z}{r}^2) \end{gathered}$$

(C15)$$=\frac{1}{P}\iint \mathrm{d}\mathit{r}\mathit{r}{|\psi (\mathit{r},0)|}^2+\frac{\lambda z}{i2\pi P}\iint \mathrm{d}\mathit{r}{\psi }^{*}(\mathit{r},0)\nabla \psi (\mathit{r},0).$$

The first integral is a centroid measurement, and from Eq. (C4) we recognize the second integral as a C tilt, giving us

(C16)$$\mathit{X}(z)=\mathit{X}(0)+z{\mathit{\theta }}_{\mathrm{C}}(0).$$

Thus, under vacuum propagation, the centroid moves by an amount given by the propagation distance times the C tilt, without our needing to know anything about how the phase or irradiance profiles of the beam actually change.

D. Integrating C Tilt and Centroid Motion Along a Path

Consider a beam that starts with no C tilt and a centroid at the origin, propagating from $z=0$ to $L$. If we consider an infinitesimal slice of atmosphere $\mathrm{d}z$ at position $z$, it will impart some small phase to the propagating beam. Given the irradiance profile of the beam as it passes through this portion of the atmosphere, we can apply the C-tilt operator to that new phase and calculate the differential C tilt that it applies to the beam, $\mathrm{d}{\mathit{\theta }}_{\mathrm{C}}$. Because propagation itself does not change the C tilt and C tilts add linearly, we merely add this new C tilt to the beam and can write

(C17)$${\mathit{\theta }}_{\mathrm{C}}(z)={\int }_0^z\mathrm{d}{z}^{\prime }\frac{\mathrm{d}{\mathit{\theta }}_{\mathrm{C}}}{\mathrm{d}z}{|}_{z=z^{\prime }}.$$

The final result of the previous section, Eq. (C16), can be written as the differential

(C18)$$\frac{\mathrm{d}\mathit{X}}{\mathrm{d}z}={\mathit{\theta }}_{\mathrm{C}}(z),$$

which we integrate along the path to give

(C19)$$\mathit{X}(L)={\int }_0^L\mathrm{d}z{\mathit{\theta }}_{\mathrm{C}}(z)={\int }_0^L\mathrm{d}z{\int }_0^z\mathrm{d}{z}^{\prime }\frac{\mathrm{d}{\mathit{\theta }}_{\mathrm{C}}}{\mathrm{d}z}{|}_{z=z^{\prime }}.$$

By the Cauchy formula for repeated integration, this becomes

(C20)$$\mathit{X}(L)={\int }_0^L\mathrm{d}z(L-z)\frac{\mathrm{d}{\mathit{\theta }}_{\mathrm{C}}}{\mathrm{d}z}.$$

This result shows that the slice of atmosphere at each position along the path moves the centroid at the end of the path by the C tilt applied by that slice times the lever arm over which the tilt acts.

Funding

Joint Directed Energy Transition Office (Air Force Contract FA8702-15-D-0001).

Acknowledgment

Distribution Statement A. Approved for public release. Distribution is unlimited. Any opinions, findings, conclusions, or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the government of the United States.

REFERENCES

1. G. A. Tyler, “Bandwidth considerations for tracking through turbulence,” J. Opt. Soc. Am. A 11, 358–367 (1994). [CrossRef]  

2. C. Hogge and R. Butts, “Frequency spectra for the geometric representation of wavefront distortions due to atmospheric turbulence,” IEEE Trans. Antennas Propag. 24, 144–154 (1976). [CrossRef]  

3. L. C. Andrews and R. L. Phillips, Laser Beam Propagation through Random Media, 2nd ed. (SPIE, 2005).

4. R. J. Sasiela, Electromagnetic Wave Propagation in Turbulence, 2nd ed. (SPIE, 2007).

5. S. M. Rytov, Y. A. Kravtsov, and V. I. Tatarskii, Principles of Statistical Radiophysics (Springer-Verlag, 1989), Vol. 4.

6. V. I. Klyatskin and A. I. Kon, “On the displacement of spatially-bounded light beams in a turbulent medium in the Markovian-random-process approximation,” Radiophys. Quantum Electron. 15, 1056–1061 (1972). [CrossRef]  

7. J. W. Goodman, Introduction to Fourier Optics, 4th ed. (Roberts & Company, 2017).

8. H. F. Talbot, “Facts relating to optical science, No. IV,” Philos. Mag. 9(56), 401–407 (1836). [CrossRef]  

9. L. Rayleigh, “On copying diffraction-gratings, and on some phenomena connected therewith,” Philos. Mag. 11(67), 196–205 (1881). [CrossRef]  

10. M. Mansuripur, Classical Optics and its Applications, 2nd ed. (Cambridge University, 2009).

11. B. J. Herman and L. A. Strugala, “Method for inclusion of low-frequency contributions in numerical representations of atmospheric turbulence,” Proc. SPIE 1221, 183–192 (1990). [CrossRef]  

12. J. W. Goodman, Statistical Optics (Wiley Classics, 2000).