Long-range propagation through inhomogeneous turbulent atmosphere: analysis beyond phase screens

Many applications rely on the propagation of electromagnetic waves through extended regions of the atmosphere over which the refractive index can vary in a complex manner. Gradients and curvature of the mean refractive index profile result in ray bending and the associated phenomena of mirages, atmospheric lensing, and wave trapping in parabolic cavities. Stochastic refractive index fluctuations due to turbulence cause a random displacement of the trajectory and give rise to the wander, or spot dancing, of a propagating optical beam. In this paper we model these features of the refractive index profile locally and describe propagation through the corresponding regions. We derive formulas for the mean ray path that capture the effects of both atmospheric turbulence and variations in the mean refractive index profile, including the non-paraxial effects associated with the bending of the guiding ray path. We also give formulas for the mean-squared transverse displacement of a ray from the mean trajectory, which can provide for example an estimate of the magnitude of the beam wander due to turbulence.


Introduction
Many technologies that rely on electromagnetic (EM) wave propagation through various regions of the atmosphere are significantly affected by the multiscale variations in the atmospheric medium. EM waves propagating through the inhomogeneous atmosphere are acted on by both mean-field refractive index gradients and stochastic refractive index fluctuations due to turbulence. Atmospheric turbulence deforms the wavefront as sections passing through regions where the refractive index is larger travel more slowly than those moving through regions where the refractive index is smaller. In this way turbulence locally focuses and defocuses the wavefront which results in intensity fluctuations, referred to as scintillation. Turbulence also causes an optical beam propagating through the atmosphere to wander randomly in a given observation plane transverse to the propagation direction, a phenomenon called beam wander or spot dancing. Beam wander also contributes to scintillation [1,2], and is one of the main causes of major power loss in Free Space Optical (FSO) communication links [3]. Scintillation limits the angular resolution of ground-based telescopes, causes speckles in images [4], and can increase the link error probability in optical communication systems [5]. In astronomical photometry, intensity fluctuations due to atmospheric turbulence may be larger than the size of the variations that must be detected in order to study, e.g. extrasolar planets [4]. Scintillation caused by the scattering of radio waves by plasma density irregularities in the ionosphere can significantly impact the performance of transionospheric navigation systems like the Global Positioning System (GPS) [6][7][8][9]. As well as affecting GPS microwave signals, ionospheric turbulence also causes distortions in spaceborne synthetic aperture radar (SAR) images [10]. In addition to the random fluctuations due to turbulence, mean variations in the atmospheric properties also have a substantial effect on propagation. The bending of radiowave trajectories due to the altitude dependence of the mean refractive index must be accounted for in, for example, airborne radar calculations [11].
The theoretical understanding of the effects of the turbulent atmosphere on EM wave propagation grew out of Kolmogorov's 1941 theory of turbulence [12]. Obukhov [13], Corrsin [14], Tatarskii [15], and others built on Kolmogorov's theory to develop models for temperature and refractive index fluctuations. Keller [16] was one of the first to study the classical equations of wave propagation in the situation of a random medium, and Tatarskii was instrumental in unifying the theory by combining these equations with Kolmogorov's description of the medium fluctuations due to turbulence. Various approaches, many of which can be classified as perturbation methods, have been developed for analyzing the resulting equations [15][16][17][18][19][20]. Much credit should also be given to Andrews and Phillips for their work in building off this more general framework to develop a comprehensive theory for laser beam propagation through turbulence [21].
A common approach to studying propagation through turbulence is to use the method of random phase screens. A phase screen is an artificial thin turbulent layer along a propagation path that is used to approximate a random medium that extends between a transmitter and a receiver. When the thickness of a turbulent layer is sufficiently small, only the phase, and not the amplitude, of an optical wave will be significantly affected by the layer. By placing a phase screen at a certain distance between the transmitter and receiver, values of some statistical quantities associated with the wave can be made to match those of a wave that propagates through an extended random medium [21]. Wave optics simulations often use multiple random phase screens placed along the propagation path, and these are usually equally-spaced and generated based on the Kolmogorov statistics (see, e.g. [22]). However, conventional phase screen approaches fail in critical regions because they use piecewise constant approximations for the mean refractive index profile. Such approximations do not sufficiently account for gradients of the refractive index field and also fail in the vicinity of local extrema where locally quadratic approximations are required to resolve EM wave fields. Because of this, they fail to capture important EM wave propagation phenomena, such as mirages, lensing, caustics, and trapping of EM waves in parabolic cavities, that are due to the inhomogeneity of the medium coupled with stochastic fluctuations. In addition, in both phase screen methods and extended random medium models, the structure of the stochastic fluctuations is usually modeled as homogeneous, isotropic turbulence having an equilibrium inertial range in its spectrum. These assumptions are not realistic in atmospheric applications when inhomogeneity and anisotropy are strong due to specific production mechanisms such as jet streams and gravity waves [23][24][25][26].
In this paper we present a method for describing longrange propagation through inhomogeneous, multiscale, and stochastic variations in the atmospheric refractive index. This paper extends the work of [27,28], in which we developed a novel approach consisting of the following steps: (i) identify critical regions containing large gradients, curvature, or local extrema of the refractive index profile (ii) locally approximate the refractive index using linear, quadratic, and stochastic terms These local approximations can capture, for example, scintillation-producing irregularities such as those formed at the edges of turbulent Rayleigh-Taylor ionospheric plasma bubbles and atmospheric layers [8]. This method is able to resolve the non-paraxial effects associated with the bending of the guiding ray path, including the phenomena of mirages and focusing/defocusing in long-range EM wave propagation through the inhomogeneous, turbulent atmosphere.
In [28] we considered the problem of characterizing the ray trajectory in a random medium with a mean refractive index gradient. Such a gradient results in the phenomenon of a mirage, where an object's observed location is displaced from its actual location. We derived formulas for the mean ray path in both the situation of isotropic stochastic fluctuations and an important anisotropic case. Our complementary work [27] studied propagation in the presence of both curvature and stochastic fluctuations in the refractive index profile. Curvature of the refractive index profile can cause the focusing or defocusing of light rays or a propagating beam. Such lensing can have detrimental effects on, for example, imaging through the atmosphere, where it can cause features to appear stretched or regions to appear artificially bright. We derived formulas for the mean ray path in such a region and also investigated the case of a local maximum of the mean refractive index profile, which can act as a gradient-index lens. For this situation, we showed that in the presence of stochastic fluctuations the focal length associated with the mean trajectories is smaller than in the deterministic setting.
The organization of the paper is as follows. In section 2, starting with Maxwell's equations we briefly review the mathematical description of EM wave propagation. In section 3 we discuss the nature of both the deterministic and the stochastic refractive index variations in the atmosphere. In section 4 we describe the refractive index models that we use for different regions in the atmosphere. These models account for both the variations in the mean refractive index as well as the random refractive index fluctuations due to turbulence. In section 5 we derive formulas for the mean trajectories in regions described by these refractive index models. In section 6 we derive formulas for the mean-squared magnitude of the transverse displacement of a ray from the mean trajectory. We discuss conclusions in section 7.

EM wave propagation
The electric E and magnetic H fields satisfy Maxwell's equations: The permittivity ε, the permeability μ, the charge density ρ, and the electric conductivity σ are properties of the medium. The equations (1a) and (1b) are Gauss's law and Gauss's law for magnetism, and equations (1c) and (1d) are Faraday's and Ampère's laws, respectively. In the atmosphere all of the properties of the medium that enter into Maxwell's equations can depend on both position and time. While accurate numerical methods have been developed for solving Maxwell's equations (see, e.g. [29,30]), simpler approximate (i.e. reduced) models are necessary for many applications that require fast solutions over very large regions. We assume here that the medium is time-independent, charge-free (ρ = 0), and lossless (σ = 0). We further assume that the permeability μ = μ 0 is constant. We express the permittivity as e e e = ( ) ( ) r r r 0 where ε 0 is the permittivity of free space and the dielectric constant ò r depends on position. We assume that the electric and magnetic fields are time-harmonic, that is, they can be expressed in the form Here ω is the frequency and  and  are the time-independent parts of the electric field and the magnetic field, respectively. Then Faraday's law (1c) and Ampère's law (1d) take the time-independent form Taking the curl of (2a) and then using (2b) and (1a) with ρ = 0 we obtain the following equation for  : . Using these relationships, In free space where the refractive index is constant the third term on the left-hand side of (4) is equal to zero. In an inhomogeneous medium this term is nonzero and describes the coupling between the three components of  and the resulting polarization effects (see, e.g. [31]). Dropping this term decouples the equations for the components of  and therefore leads to a substantially simpler model that can be adequate for applications in which polarization is not measured. The resulting model is the vector Helmholtz equation In considering this equation, it is common and useful to look for solutions of the form where the amplitude (or envelope)  and phase Ψ=k 0 ψ are real-valued. In many propagation scenarios of interest, especially in the case of visible light, the wavelength l p = k 2 0 is small compared to the other length-scales in the problem. In this situation the wavenumber is large and so it is useful to consider the leading-order terms in k 0 in the equation obtained after substituting (7) into (6). After matching terms of the same order in k 0 , the leading ( ) O k 0 2 equation involves only the function ψ, called the eikonal, and not the amplitude . It is therefore called the eikonal equation: Finding the field u by first solving (8) and then solving the next order equation, which depends on ψ, for the amplitude  is called the geometrical optics approximation. The characteristics of the first-order equation (8) are called rays, as they correspond to the notion of light rays. From (8) one can see that a ray is orthogonal to the level surfaces of the solution of the eikonal equation. That is, under the geometrical optics approximation a ray is perpendicular to the wavefronts of the solution of the Helmholtz equation (6). Letting s denote arclength, a ray = ( ) ( ( ) ( ) ( )) r s x s y s z s , , satisfies the ray equation Let y g.o. denote the solution of (8) so that y k 0 g.o. is the geometrical optics approximation of the phase of the solution of the Helmholtz equation (6). The difference between the values of y g.o. at two points r A and r B can be found by integrating the refractive index along the connecting ray (see, e.g. [32]): Here r is the ray connecting r A and r s , B is arclength along this ray, = ( ) r r s A A , and = ( ) r r s B B . Any description of EM wave propagation requires information about the properties of the medium. The complete description given by Maxwell's equations requires the permittivity, permeability, charge density, and electric conductivity to be known. On the other hand, assuming the wavenumber k 0 is known, the solution of the Helmholtz equation (6) and that of the ray equation (9) only depend on the medium through the refractive index. Therefore, information about the refractive index alone can give a description of propagation that is sufficient for many applications. We turn now to the problem of modeling the refractive index in the atmosphere.

Refractive index variations in the atmosphere
Variations in the refractive index in the atmosphere are caused by inhomogeneities of a wide range of sizes. These include both variations in the mean refractive index and refractive index fluctuations due to turbulence that cannot be accounted for by deterministic models. Larger-scale mean variations in atmospheric properties, such as the decreases in pressure, temperature, and humidity over altitude increases on the order of kilometers, can often be resolved. It is an oversimplification, however, to assume that changes in the mean refractive index always correspond to larger-scale variations. For example, there is evidence that curvature features of certain sections of the vertical refractive index profile on the order of meters occur consistently at the same time of day and therefore may be considered deterministic and predictable [33,34]. On the other hand, turbulence gives rise to refractive index fluctuations over a wide range of scales. In the troposphere, tropopause, and lower stratosphere, turbulence mixes layers of different temperatures and densities [23,25,26] and results in inhomogeneities that range in size from about 1cm to several kilometers [35]. Also in these regions, very strong temperature gradients can occur within 'sheets' with thicknesses on the order of meters [36]. In the ionosphere, plasma instabilities result in turbulent structures of scales that range from larger than 10 5 m to smaller than 0.1 m [8,37]. In this region, the refractive index depends on the electron density [10], whose order of magnitude varies with altitude throughout the ionosphere [37][38][39], and which varies with latitude and longitude as well [7,9].
In the troposphere, variations of the mean refractive index in the vertical direction are much larger than those in the horizontal directions. Therefore, models for the troposphere commonly take the mean refractive index to depend only on altitude. In the atmosphere the refractive index differs from a value of one by a very small amount. Because of this, changes in the refractive index are commonly expressed in terms of the refractivity N which is defined as = -( ) N n 10 1 6 . Over vertical distances on the order of kilometers in the troposphere, average radiowave refractivity is often modeled as an exponentially decreasing function of altitude due to the changes in pressure, temperature, and humidity [11,40]. That is, average radiowave refractivity is modeled by where N s is the refractivity at the Earth's surface, x is height in kilometers above the Earth's surface, and H is the scale height. Typical midlatitude values for the parameters N s and H are N s =315 and H=7.35 km [11]. At low heights, e.g. less than a kilometer, a linear approximation is often used for the vertical refractivity profile [40]. At midlatitudes a representative average vertical gradient of the radiowave refractivity in the first kilometer is −40 km 1 [40]. Using this value gives the linear approximation for the radio refractive index, where n s is the value of the refractive index at the Earth's surface. Rays that propagate in such a region with a refractive index gradient are bent toward the direction of the gradient and therefore travel in curved paths. This leads to a vertical translation of a scene in images as well as the associated phenomenon of a mirage where an object's observed location is different from its actual location. When light reflected by an object reaches an observer via such a curved path, the apparent direction of the object is the direction from which the light ray enters the observer's eye. Thus, the object is observed to be located at a position that is displaced in the direction of the negative gradient of the refractive index from the object's actual location. For example, in a typical setting, the mean vertical refractive index gradient near the ground results in a vertical displacement of several meters when an object is observed at a horizontal distance of around 13-15 km [41]. A linear model such as (11) of the refractive index is sufficient to capture this ray bending and the resulting phenomenon of mirages.
A linear model cannot, however, describe the local curvature of the refractive index profile, nor can a monotonically decreasing model resolve local extrema. A local maximum of the profile can act as a gradient-index lens that focuses a beam propagating around the height of this maximum. Conversely, a local minimum defocuses propagating rays and thereby spreads a beam. Such focusing and defocusing is sometimes referred to as atmospheric lensing. This can significantly affect imaging applications in particular because focusing and defocusing can cause a stretching and a compression, respectively, of a scene in images. When measurements are used to generate plots of the refractive index in the troposphere as a function of altitude on a 1 km scale, the profile may exhibit multiple local minima and maxima [42,43]. A local maximum can arise, for example, near a weather front, where it results from the slope of the front and the presence of smaller vertical refractive index gradients on the cold side than on the warm side [44]. More generally, atmospheric lensing can result from not only local extrema of the refractive index profile but also any segment of the profile that has nonzero curvature. As mentioned above, one may be able to predict and resolve the curvature of certain parts of the refractive index profile. In [33,34], the value of the curvature of the vertical refractive index profile in a region within about ten meters of the ground was deduced from measurements of the resulting apparant stretching of a building in time-lapse photographic images. While the profile changed significantly over the course of a day, the largest stretch for each day typically occurred at the same time of day, and the corresponding value of the curvature parameter was consistent over a five day period. The curvature of this section of the profile resulted in an apparant stretch of about 5 m of a 33 m tall building in images taken at a distance of about 15 km.
Any theoretical description of the refractive index fluctuations due to atmospheric turbulence must be statistical in nature. The spatial structure of these fluctuations is often described in terms of either the structure function D n or the covariance function B n of the refractive index. These are defined as where E denotes expectation. The spatial structure is also commonly expressed in terms of the spectrum of the fluctuations, which is related to the covariance function through the Wiener-Khinchin theorem when the fluctuations are statistically homogeneous. The classical Kolmogorov theory of turbulence assumes that turbulence is both statistically homogeneous and isotropic [45]. Under these assumptions the values of the structure function and covariance function at ¢ ( ) r r , depend only on the distance = -¢ | | r r r between the two points, and furthermore they are related by The Kolmogorov theory gives an explicit form for the structure function and thereby also the covariance function of the refractive index. It says that there exists an inertial subrange of length-scales of the fluctuations, bounded above by the outer scale L 0 and bounded below by the inner scale ℓ 0 . The outer scale is often identified with the correlation length of the refractive index fluctuations due to turbulence. The lower bound ℓ 0 of the inertial subrange depends on the energy dissipation rate (per unit mass) ε and the kinematic viscosity ν. The length scale defined by these two parameters is referred to as the Kolmogorov microscale h K , which is given by h n e = ( ) K 3 1 4 . In the atmosphere the relationship between the inner scale and the Kolmogorov microscale is typically taken to be ℓ 0 =7.4η K [21,46]. Kolmogorov's theory states that the structure function is proportional to r 2/3 in the inertial subrange, which leads to the following form for the covariance function of the refractive index (see, e.g. [21,35]): where s n 2 is the variance of the refractive index. The two-thirds power law for the structure function is equivalent to the statement that the three-dimensional spectrum Φ(κ) is proportional to κ −11/3 over the inertial subrange. This corresponds to a one-dimensional spectrum with a −5/3 power law, which is called the Kolmogorov spectrum. Non-Kolmogorov spectra, i.e. those with a power law different from −5/3, have also been considered, see, e.g. [47][48][49][50].
There has been much work on measuring the values of the turbulence parameters C ℓ , , , where the largest values are seen in turbulent layers located at the top boundaries of clouds. Several models of the altitude dependence of C n 2 have been developed based on experimental data. One of the most widely used is the Hufnagel-Valley model, see [52,53]. Turning now to the inner and outer scales, the inner scale ℓ 0 is typically on the order of 1-10 mm near the ground and around 1 cm or larger in the troposphere and stratosphere [4,21,35]. Within the atmospheric boundary layer the outer scale is commonly taken as L 0 =0.4h where h is height above the ground. In the free atmosphere the outer scale L 0 is about 100 m in the vertical direction and a few kilometers in the horizontal direction [35]. Thus, the correlation length of the refractive index fluctuations due to atmospheric turbulence is highly anisotropic. One of the main causes of this is gravity waves, which generate strongly anisotropic refractive index fluctuations with horizontal scales of many kilometers [54].

Refractive index model
In view of the discussion of the previous section, both the qualitative nature of the mean refractive index profile and the statistical properties of the random refractive index fluctuations due to turbulence can vary substantially over the large regions of the atmosphere through which EM signals propagate. In attempting to describe the refractive index over the entire propagation region it is therefore useful to model the refractive index locally in sections over which the important features of the mean profile can be resolved and the statistical properties of the random fluctuations are approximately constant. A model for the entire region can then be constructed as a composition of these local models. In this section we describe local models for the qualitatively different sections of the refractive index profile.
We will from here on focus on the scenario of nearlyhorizontal propagation in the troposphere, tropopause, and stratosphere, which is important for many applications. For example, in airborne SAR a scene on the ground is imaged by emitting successive radiowave pulses from an antenna mounted on an airplane. The airplane may typically fly at an altitude of 5 km and be located at a distance of around 100 km from the scene being imaged [55]. In this situation, the grazing angle between the antenna and the target scene is around 3°so that each signal propagates nearly horizontally yet at the same time travels a total vertical distance over which the vertical refractive index profile may exhibit many different features. Since the mean refractive index mainly depends on altitude, each signal propagates in a direction approximately perpendicular to the direction of the most significant variations in the mean refractive index. It is in this situation that the variations in the mean refractive index have the largest effects on the propagation trajectory.
In regions in which the curvature of the refractive index profile is negligible, the mean profile can be described locally by a linear term in the refractive index model. To account for refractive index fluctuations due to turbulence it is necessary to include a stochastic term in the model. Therefore, in these regions we use the following model n ℓ for the refractive index: Here n 0 is a reference value of the refractive index which is very close to 1 in the atmosphere and β ℓ is a dimensional parameter. We assume 0ò ℓ = 1 so that ò ℓ is a small dimensionless parameter that represents the fact that the refractive index variations in the atmosphere are very small compared to the reference value. The stochastic fluctuations are represented by the mean-zero random field η which we describe in detail below. Such a model was used in [28] to describe the effects on the propagation trajectory of a mean gradient and stochastic fluctuations of the refractive index. As an example of the size of the coefficient ò ℓ β ℓ , using the representative value given in (11)  In regions where the curvature of the profile is significant, or near a local extremum of the mean profile, a different model should be used. Near an extremum, the mean refractive index function can be Taylor-expanded locally about the critical point. In this expansion the linear term is equal to zero and the quadratic term is sufficient to capture the lensing effect. Furthermore, a quadratic term in the model can be used to describe the trajectory locally near a curved segment of the mean refractive index profile. The coefficient of the quadratic term can be tuned as a parameter to match the curvature, while the rayʼs initial vertical position in the equation (9) can be used as another parameter to match the first derivative of the true profile near the ray. Therefore in these regions we use the following model: Here, ò q is a small dimensionless parameter and β q is a dimensional parameter. A negative value of β q represents a segment of the mean refractive index profile with negative curvature, e.g. a local maximum, which will focus a propagating beam. Conversely, β q >0 corresponds to a segment with positive curvature that gives rise to negative lensing (defocusing). A model of this form was used in [27] to study lensing in the presence of stochastic fluctuations. A quadratic term in the refractive index model has also been considered in the paraxial wave equation [56]. As an example of the size of the coefficient ò q β q we refer to the results of [33,34] where the value of the curvature of the vertical refractive index profile in a region within about ten meters of the ground was measured. This measured value corresponds to a value of  b » -´--1.5 10 m q q 9 2 in (14). We now describe the random field η in (13) and (14) that we use to model the stochastic refractive index fluctuations due to turbulence. To model the highly anisotropic nature of the correlation length of the random refractive index fluctuations in the free atmosphere, we take the random field η to be composed of both an isotropic part and an anisotropic part. We take the anisotropic part to consist of fluctuations only in the vertical direction, that is, to correspond to random layers. This represents the fact that in the free atmosphere the correlation length of the refractive index fluctuations is much longer in the horizontal directions than in the vertical direction. Therefore we define the random field η in (13) and (14) as where η I is a homogeneous and isotropic mean-zero random field, η A is a stationary mean-zero stochastic process, and κ I and κ A are parameters that control the relative strengths of the isotropic and anisotropic random fluctuations. We will assume that η I and η A are independent (similar assumptions have been made in [54,57]). We take η I and η A to each have variance equal to one, so that the variance of the refractive index as given by (13) and (14) is where i=ℓand i=q respectively. Since η I is homogeneous and isotropic and η A is stationary, the autocovariance functions of η I and η A are of the form for all d>0. We assume that R I is differentiable and R A is twice differentiable, and we also make the physically realistic assumption Note that ¢ = ( ) R 0 0 I is satisfied by the Kolmogorov covariance function (12). By the assumed independence of η I and η A the autocovariance function of η is

Mean ray path; mirage and lensing effects
We now turn to the problem of finding the trajectory of a signal propagating through a region described by either of the refractive index models (13) and (14). When one of these models is used for the refractive index, the solution of the ray equation (9) is a random trajectory. The best estimate of the actual trajectory is therefore given by the expected or mean ray path. We will find this by first deriving an expression for the ray path corresponding to an arbitrary realization of the refractive index, and then calculating the expectation of this expression. In order to determine a specific ray's trajectory using (9), the ray's initial position and direction must be specified. As we previously noted, in the case of the model (14) with a quadratic term, the ray's initial x-coordinate can be used as a parameter to control the first derivative of the mean refractive index profile near the ray. Let b denote the ray's initial x-coordinate. Then together with the ray equation (9) we consider the following initial conditions: 1, 0, 0 1 is the unit vector in the positive x-direction. We note that = | | u 1 since s is arclength. Let θ denote the angle between u and the yz-plane. We focus on the case of nearly-horizontal propagation, which corresponds to rays for which θ≈0. In the calculations that follow, however, it is convenient to initially treat the angle θ as an arbitrary parameter, and in doing so we use the convention that θ is positive when > · u e 0 1 . Because in (13) and (14) the dimensionless parameters ò ℓ and ò q are small compared to the reference value of one, it is natural to use perturbation methods to find a ray's trajectory in regions described by these models. We use such a method to derive the mean ray path to order  i 2 in these regions. We follow the method of Keller, who derived an expression for the mean ray path in the case in which the mean refractive index is constant and the random fluctuations are isotropic [16]. We first consider the case in which the refractive index n in the ray equation (9) is defined by (13). By viewing r as a function of both s and ò ℓ , we Taylor-expand r in the variable ò ℓ about ò ℓ =0. In this way we can express ( ) r s as the following expansion in powers of ò ℓ : We substitute (18) into (9). We equate terms of the same order in ò ℓ , thereby obtaining a sequence of equations to be solved recursively for r r r , , ,... Next, we solve in succession the order-ò ℓ equation for r 1 in terms of r 0 and the order-ò ℓ 2 equation for r 2 in terms of r 1 and r 0 . We then calculate the expectations of ( ) r s 1 and ( ) r s 2 in order to find the expectation of ( ) r s to order  ℓ 2 . First, by defining  =  - ( · ) u u T so that ∇ T denotes the transverse gradient to the initial direction u, we can express r 1 as . We see that the stochastic term in the refractive index does not contribute to the expectation of r 1 . It does, however, give rise to an order- ℓ 2 contribution to the mean trajectory that appears in the expression for the expectation of r 2 . The expression for ( ) r s 2 is lengthy and so we omit it, but we describe how to calculate its expectation. The contribution to [ ( )] r E s 2 due to the stochastic term involves the quantities I T I 0 2  I T I 0 1 and for t 1 t 2 s, as well as similar expressions. We show how to calculate (22a); (22b) can be found in a similar way. To calculate (22a) we define the auxiliary functions n h , and we let  ( ) v T and  ( ) w T denote the transverse gradients to u with respect to the variables v and w, respectively. Then the expression (22a) can be written as For an arbitrary initial angle θ, the expression for the expectation of ( ) r s 2 is lengthy. However, for the case θ=0 which describes well nearly-horizontal propagation, the expression simplifies considerably. Using (19), (21), and (23), we find that for θ=0, to order  ℓ 2 , the expected position of the endpoint of a ray of length s is given by On the right-hand side of (24) the first two terms consist of the order-one straight-line trajectory r 0 . The third term is due to the linear term in the refractive index and describes the bending of the ray toward the direction of the mean gradient that gives rise to the mirage effect. The fourth term, which is also due to this linear term, is a correction to the component in the initial direction. The last two terms describe, to order  ℓ 2 , the effect of the stochastic fluctuations on the mean ray path. Since  ( ) R 0 0 A , and  ¢( ) R r 0 I as well for physically realistic models η, these terms reduce the component of [ ( )] r E s in the initial direction u. While the stochastic fluctuations contribute to the components in the directions orthogonal to u of a given realization of ( ) r s (this is seen in the order- ℓ correction r 1 given by (20)), these contributions average to zero to order  ℓ 2 when the expectation is taken. For the case in which the mean refractive index is constant (β ℓ =0) and the random fluctuations are isotropic (κ A =0), the formula (24) reduces to the one derived by Keller [16].
We now derive the mean ray path for regions described by (14). We again restrict our attention to the case θ=0 which describes well nearly-horizontal propagation. By a similar calculation as before, we find that for θ=0 the mean ray path to order ò q 2 is given by As in (24) the first two terms consist of the order-one straight-line trajectory r 0 . The third term describes the lensing effect that results from the quadratic term in the refractive index. For b < 0 q , this term causes a ray to bend in the x-direction toward the plane x=0 by an amount that depends on the ray's initial distance | | b from this plane. This gives rise to the focusing of rays. For β q >0, rays are bent away from the plane x=0, resulting in negative lensing. The next two terms are higher-order corrections due to the quadratic term in the refractive index. The last two terms, which describe the effect of the stochastic fluctuations, are the same as those in (24) (with ò ℓ replaced by ò q ).
In the case in which the random fluctuations of the refractive index are isotropic (i.e. κ A =0), the contribution of these fluctuations to the mean ray path can be expressed in terms of the structure constant, inner scale, and outer scale by using the Kolmogorov covariance function (12). For the models (13) and (14), and the case κ A =0, the covariance function B n of the refractive index is given by where i=ℓand i=q, respectively. Using (12) and (26)  For  s L 0 , the first term inside the parentheses on the righthand side dominates, so in such a case the value of L 0 does not need to be specified in order to estimate the magnitude of (27). From (27) we see that the smaller the inner scale ℓ 0 is, the larger the reduction of the component of [ ( )] r E s in the initial direction u. This is because the more rapid the medium fluctuations are, the more quickly the direction of a given realization of the ray trajectory varies.

Transverse displacement of ray due to turbulence
In the formulas (24) and (25) for the mean ray path, the stochastic fluctuations only contribute to the component of [ ( )] r E s in the initial direction u. As previously mentioned, this is because the contributions of these fluctuations to the components of ( ) r s in the directions orthogonal to u average to zero (to second order) due to cancellations when the expectation is taken. However, for a given realization of the refractive index, the fluctuations cause a displacement of the endpoint of the ray from the initial direction u (this is in addition to that caused by the linear and quadratic terms). For both of the refractive index models (13) and (14), in the case θ=0 the displacement of the endpoint of the ray from its mean value is to order-ò i orthogonal to the initial direction u (for the case κ A =0 of isotropic fluctuations, this is true for any initial angle θ). For the refractive index model (13) this can be seen from the expression (20) for the order-ò ℓ correction r 1 to the ray path.
In this section we derive formulas for the mean-squared magnitude of this transverse displacement of a ray from the mean trajectory. This can provide an estimate of the wander of a propagating optical beam due to turbulence. We will consider separately the cases of isotropic and vertical random fluctuations in the refractive index.

Transverse displacement due to isotropic fluctuations
We first consider the case in which the random fluctuations of the refractive index are isotropic, i.e. κ A =0 in (15). The component transverse to the initial direction u of the displacement of the endpoint of a ray from its mean value is given by Here I is the identity matrix and u T is the transpose of u. The expression (28) is a mean-zero random vector. To give a measure of its typical size, we calculate the second moment of its magnitude: To order- i 2 and for   = ℓ q , (29) is the same for both of the refractive index models (13) and (14). We find (29) to order- i 2 by a calculation similar to the ones by which we obtained the mean ray paths (24) and (25). We get where i=ℓ and i=q when the refractive index is given by (13) and (14), respectively. As before, we can use the Kolmogorov covariance function (12) to express (30) in terms of the structure constant, inner scale, and outer scale. We find, for s>L 0 , -- This formula expresses the mean-squared magnitude of this transverse displacement in terms of the propagation distance, structure constant, inner scale, and the dimensionless ratio ℓ 0 /L 0 which characterizes the width of the inertial subrange. For , the first term on the right-hand side dominates, giving the magnitude of (32). As an example of the magnitude of the transverse displacement due to isotropic fluctuations in the atmosphere, for C n 2 on the order of 10 −13m 2 3 , a propagation distance of around 10 km, and a typical value for ℓ 0 between a few millimeters and 1 cm, the square root of the expression given in (32) is on the order of 1 m.

Vertical displacement due to vertical fluctuations
Now we consider the case in which the random fluctuations of the refractive index are only in the vertical direction, i.e. κ I =0 in (15). We are interested in the vertical component of the displacement of the endpoint of a ray from its mean value: - Like (29), to order- i 2 and for ò ℓ =ò q , (34) is the same for both of the refractive index models (13) and (14). We restrict our attention to the case θ=0 in which the ray initially points in a horizontal direction. In this case and for this type of fluctuations, the component transverse to the initial direction u of the displacement of the endpoint of the ray from its mean value is, of course, entirely in the vertical direction, i.e. it is given by (33). By a similar calculation as before, we find, to order- i 2 , In considering the size of (35), the magnitude of ( ) R 0 A can be estimated from the value of the inner scale of the anisotropic component of the fluctuations. Typical values of this parameter in the atmosphere have been obtained from measurements, see, e.g. [54].

Discussion
We have considered the problem of describing long-range EM wave propagation through the atmosphere where both variations in the mean refractive index and stochastic refractive index fluctuations due to turbulence must be taken into account. Over extended regions of the atmosphere both the qualitative nature of the mean refractive index profile and the statistical properties of the fluctuations due to turbulence can vary substantially. While the mean refractive index in the troposphere is commonly modeled as a function of altitude only, the actual vertical profile can be complex. The mean vertical gradient varies with altitude and the profile may exhibit sections with significant curvature as well as multiple local extrema. Accounting for the gradient of the mean profile is necessary to describe the bending of the trajectory and the associated phenomenon of mirages. Furthermore, the local curvature and extrema of the mean profile must be modeled in order to capture the phenomenon of atmospheric lensing. While the random fluctuations due to turbulence are often described by the isotropic Kolmogorov model (12), the strength of these fluctuations in general varies dramatically over extended regions and the correlation length of these fluctuations is highly anisotropic.
To capture these features, we proposed local models for the qualitatively different sections of the refractive index profile. In regions where the curvature of the mean profile is negligible, the model (13) containing both a linear term and a stochastic term is used. We derived the formula (24) for the mean ray path in these regions, which describes both the bending of rays toward the direction of the mean gradient and the effect of turbulence on the mean trajectory. In regions where the curvature of the mean profile is significant, or near a local extremum, the model (14) containing a quadratic term can be used. For these sections of the profile the formula (25) for the mean ray path captures the associated phenomenon of atmospheric lensing. A guiding trajectory describing propagation over an extended region can then be found by using the formulas (24) and (25) to describe the path locally as it traverses the qualitatively different sections of the refractive index profile.
As described in formulas (24) and (25), the random refractive index fluctuations due to turbulence reduce the component of the mean trajectory in the initial propagation direction. Turbulence also causes the endpoint of a ray to be displaced in the transverse directions from the mean trajectory, which is associated with the wander of a propagating optical beam. To describe this effect, we derived the formulas (30) and (35) for the mean-squared magnitude of the transverse displacement that is due to the different types of fluctuations. The formulas (24), (25), (30), and (35) describe in a robust way the effects of the random fluctuations on the trajectory since we have only made very general assumptions about the random field that models these fluctuations. For the case in which the random refractive index fluctuations due to turbulence are statistically isotropic we have used the Kolmogorov covariance function to express the effects of the fluctuations on the trajectory in terms of the standard turbulence parameters. This enables estimates of the magnitudes of these effects in different atmospheric propagation scenarios.