Non-Classic Atmospheric Optical Turbulence: Review

Theoretical models and results of experimental campaigns relating to non-classic regimes occurring in atmospheric optical turbulence are overviewed. Non-classic turbulence may manifest itself through such phenomena as a varying power law of the refractive-index power spectrum, anisotropy, the presence of constant-temperature gradients and coherent structures. A brief historical introduction to the theories of optical turbulence, both classic and non-classic, is first presented. The effects of non-classic atmospheric turbulence on propagating light beams are then discussed, followed by the summary of results on measuring the non-classic turbulence, on its computer and in-lab simulations and its controlled synthesis. The general theory based on the extended Huygens–Fresnel method, capable of quantifying various effects of non-classic turbulence on propagating optical fields, including the increased light diffraction, beam profile deformations, etc., is then outlined. The review concludes by a summary of optical engineering applications that can be influenced by atmospheric non-classic turbulence, e.g., remote sensing, imaging and wireless optical communication systems. The review makes an accent on the results developed by the authors for the recent AFOSR MURI project on deep turbulence.


Introduction
The aim of the introductory section is three-fold: first, outlining the historical development of concepts and parameters characterizing 3D turbulent motion; second, thoroughly familiarizing the reader with the classic theory of turbulence; and third, briefly discussing the 2D turbulence theory being a stand-alone, comprehensive example of non-classic turbulence.

Early Studies
Leonardo da Vinci was one of the first scholars who attempted to visualize, classify and comprehend the phenomenon of turbulence. Figure 1 shows the drawing that he made in 1509 on recording the experimental observation of water surface turbulence produced by differently shaped objects. The appearance of the turbulent eddies behind the flat obstacle is evident. It is currently argued that Leonardo could have been exposed to ancient books that flooded Italy after the fall of Constantinople earlier in 1453. This suggests that the phenomenon of turbulence could have been appreciated for millennia.
In the middle of the 18th century, Leonard Euler established the broad field of fluid dynamics by introducing a set of partial differential equations governing the motion of incompressible, adiabatic and inviscid flows. The Euler equations are the consequence of the balance between the fluid's energy and momentum and incorporate the continuity law (conservation of mass) [1]. Almost a century later, Gabriel Stokes derived a set of equations that generalized Euler's equations to viscous fluids and are now known as the Navier-Stokes equations [2]. By the middle of the 19th century, the interest in turbulence resurfaced. On studying instabilities in the pendulum's motion, G. Stokes introduced a parameter that characterizes fluid flows by the amount of their mixing [3], currently known as the Reynolds number (after Osborn Reynolds who was the first to make practical use of the concept [4]). The Reynolds number characterizes the ratio of inertial to viscous forces in the fluid: where u [m/s] is the speed, L [m] is the characteristic linear dimension, and ν [m 2 /s] is the kinematic viscosity of the fluid. The transition from laminar (mixing-free) to turbulent flows depends on geometry-from 2.3 × 10 3 to 5.0 × 10 4 for pipe flow to 10 6 for boundary layers. In atmospheric cyclones, the number can reach values on the order of 10 12 .
On developing the tools for weather forecasting, Lewis Fry Richardson came up in the early 1920s with the construction that is now known as the direct energy cascade or Richardson energy cascade [5]. In a classic 3D turbulent flow, the energy is mechanically injected into a volume filled with a fluid by means of a velocity field, such as wind in a gas or a current in a liquid. If the initial distribution of some scalar field (temperature, pressure, concentration of a chemical compound) within the volume is not uniform, the injected energy creates a few sizable eddies, which later break down into a larger number of smaller ones. The process repeats until a sufficiently small scale is reached and the mechanical energy starts dissipating into heat. A turbulent eddy is intuitively understood as a region of space in which the scalar field remains fairly correlated. In 1925, Ludwig Prandtl introduced a parameter, currently known as the Prandtl number, which characterizes the ratio of momentum diffusivity to diffusivity of an advected scalar. It is defined by the following expression: [6] Pr = ν/α j , where α j [m/s 2 ] is the diffusivity relating to the nature of the advected scalar: thermal for temperature advection, mass for molecular advection, etc. This parameter plays a crucial part in characterizing non-classic turbulent regimes. For atmospheric turbulence, the value of Prandtl number is about 0.72. Having been greatly influenced by the work of Richardson, in 1935, Sir Geoffrey Ingram Taylor introduced the concept of locally isotropic 3D turbulence [7], the isotropy (or invariance, in the limit) being understood as the weak dependence of the flow's characteristics on orientations in space. Taylor is also known for his frozen turbulence hypothesis, suggesting that all the eddies are advected by the mean streamwise velocity without changes in their properties. Using the conception of Taylor, Theodore de Karman and Leslie Howarth introduced the 3D correlation tensor of the turbulent flow and reduced it to a scalar correlation under the isotropy constraint [8]. Theory (1941Theory ( -1942 After establishing the rigorous, mathematical theory of random functions in 1930s [9,10] Andrey Nikolaevich Kolmogorov and his students attempted its application to isotropic turbulence. Greatly influenced by the results of de Karman and Howarth, Millionshchikov applied the concept of Kolmogorov's ensemble averages to a homogeneous and isotropic turbulent flow [11]. He also developed the system of equations governing the velocity correlation and implemented quasi-normal approximation to obtain the closed system of the second-and third-order correlations [12].

Kolmogorov
Between 1941 and 1942, Kolmogorov established the foundations of the classic turbulence theory by formulating three hypotheses regarding the equilibrium regime in a homogeneous turbulent flow at sufficiently high Reynolds numbers and for sufficiently small scales [13][14][15][16][17] (see also important parallel work by Obukhov [18][19][20]). This theory largely elaborated on the Richardson turbulent energy cascade by assigning to a turbulent eddy a scale l, a velocity vector u(l) and a time τ(l). It was then conjectured that the cascade is formed by the largest eddies, with a scale L 0 (integral or outer scale), L 0 ≤ L, L being the scale of entire region of mixing. It can be estimated as follows: where K is the turbulent kinetic energy (reduced to unit mass) as follows: and ε [m 2 s −3 ] is the energy dissipation rate. Then, by defining the turbulent Reynolds number as and applying it to the integral scale, one obtains the following: attaining very large values. The first Kolmogorov hypothesis states: (i) the turbulent motion of eddies with scales l << L 0 is statistically isotropic. It appears possible to estimate a threshold scale L I that separates the regimes with isotropic (L I < l < L 0 ) and anisotropic (l < L I ) eddy scales, as L I ∼ L 0 /6. The second hypothesis postulates that (ii) the statistics of the small-scale (l < L I ) eddy motion have a universal form that is uniquely determined by ε and ν. This range of scales is also known as the universal equilibrium range. It follows from (ii) that a set of parameters solely depending on ε and ν can be introduced: These are the characteristics of the smallest eddies before dissipation, carrying Kolmogorov's name. Note that Rt starts at a very large number when energy is injected at scale L 0 and decreases to unity at scale l K . The third hypothesis states that (iii) the statistics of motion of eddies with scales in the range l K << l << L 0 have a universal form that is uniquely determined by ε and are independent of ν.
There is another scale, l D ≈ 60l K , that splits the universal dissipation range (l < l D ) in which motion is the result of viscous forces.
The Kolmogorov theory established the regime of classic turbulence conceptually, and also led to specific mathematical law governing the structure of turbulent velocity correlation function or, equivalently, its energy spectrum. On introducing spatial frequency as κ = 2π/l (or κl K for dimensionless version) and energy spectrum E K (κ) such that one finds at once that, in order for E K (κ) to depend on κ and ε, i.e., to be in the form it suffices to perform dimensional analysis.
= m 3 s −2 into the equation above, one finds that p = 2/3 and q = −5/3, thereby leading to the famous Kolmogorov law as follows: where C = 1.5 is a data-fitted constant. One may also define the dissipation rate spectrum, D K (κ) [m 3 /s 3 ], as follows: It was estimated that most of the energy (80%) is contained in the energy production range L I < l < L 0 , and most of the dissipation (90%) occurs in the interval 8l K < l < l D . Additionally, in terms of a temporal scale of an eddy, on converting to its lifespan τ = K/ε, 90% is spent in the production range and 10% in the inertial range.

Batchelor-Leigh-Kraichnan Theory
In the late 1960s, Robert Harry Kraichnan [21], Cecil Leith [22] and George Batchelor [23] (see also earlier important work [24]) developed the theory of 2D turbulence. Unlike in the 3D case, in addition to energy K, the turbulent flow must also be characterized by enstrophy (reduced to unit mass): where w = − → ∇ × u is the flow helicity. The power spectra of energy E K and enstrophy E Z can be shown to be in a simple relation as follows: Unlike in 3D turbulence, described by the direct turbulent energy cascade, in which the energy is injected at the largest scale and dissipates at the lowest scale, the 2D turbulence is characterized by a dual cascade: the enstrophy undergoes the direct cascade from the intermediate injection scale to the lowest scale, having the κ −3 power law and the energy exhibits the inverse energy cascade from the intermediate scale to the largest scale, having the same (in form) power law κ −5/3 as the classic 3D turbulence.
Even though the 2D turbulence model could not be directly applied to 3D turbulence research and applications, it provided the invaluable insight into both classic and non-classic turbulent regimes and helped in describing various equilibrium and nonequilibrium turbulent regimes, particularly those appearing in rotating and conducting fluids. We refer the reader to a recent comprehensive review of these studies [25].

Measurements and Synthesis
In this section, we bring together the results of some early and recent measurements, both thermodynamic and optical, as well of computer simulations, indicating the possibility of air turbulence to deviate from its classic regime, as originally postulated by Kolmogorov.

Anisotropic Turbulence
The classic Kolmogorov theory, based on Richardson's energy cascade, implies that in the inertial sub-range of scales, all turbulence eddies are isotropic. However, it appears possible, in free atmosphere, within the stably-layered stratosphere, for the turbulence to become anisotropic (direction-dependent) at large spatial scales (e.g., [26,27]). A similar situation may appear in the close proximity of a hard boundary, i.e., ground surface, a building, an airplane, etc., i.e., within a boundary layer [28][29][30][31].
In their 1970 experimental paper, Consortini, Ronchi and Stefanutti illustrated in the laboratory (over 130 m path, at 1 m height) that in the close proximity of a boundary (ground), a pair of laser beams, co-aligned and set to propagate along the horizontal direction, would lead to short-exposure relative beam wander ("relative dancing"), depending on the placement of the two sources [32]. In particular, the statistics were shown to qualitatively and quantitatively differ for horizontal and vertical orientations, implying the presence of turbulent anisotropy. This measurement procedure suggested a very simple, optics-based method for imprinting the signature of air turbulence onto the intensity statistics of propagating light. Since the fluctuating temperature of the turbulent air is the dominating random process affecting the fluctuating refractive index, the air velocity turbulence can be directly related to the refractive index turbulence, also termed optical turbulence. In particular, because temperature fluctuations are largely responsible for optical turbulence of the air, the anisotropic behavior of the temperature field carries over that in the fluctuating refractive index field.
In order to develop the theory, a 3D Gaussian (ellipsoidal) correlation function model of the refractive index fluctuations of the turbulent air was implemented: where n 2 is the mean-square value of the refractive index about its mean value, and ∆j 0 is a typical deviation in direction j, (j = x, y, z). Physically, the anisotropic turbulent medium can be viewed as the random collection of ellipsoidal lenses with average semiaxes lengths and orientations corresponding to geometrical directions x, y (horizontal) and z (vertical). In the vicinity of a horizontal boundary, for example, ground, ∆z 0 < ∆x 0 , ∆y 0 . This discrepancy leads to qualitatively different refraction scenarios along the horizontal and vertical axes. The dancing of the two beams is characterized by the wave structure function, being a variance of the difference of the optical field measured at several separation distances. The structure functions D y and D z , measured for horizontal and vertical sources' placement, are shown in Figure 2A,B, respectively. The solid lines denote the theoretical joint structure functions, while dots denote the experimental data. The discrepancy in D y and D z is related to the shape of B n along directions x and y.
Since the 1990s, high-altitude atmospheric temperature field measurements have become available. A strong anisotropy in temperature was first found by Grechko et al. [33] from experimental observations of starlight scintillation at the intermediate altitudes.
Another experimental campaign by Dalaudier et al. [34] revealed the presence in the atmosphere, at altitudes of up to 25 km, of very strong positive temperature gradients within very thin layers, or sheets, with vertical distortions up to 10 m and of horizontal extensions larger than 100 m. The anisotropy in the experimentally measured optical wavefront tilt's statistics was also revealed by Belenkii: the outer scale in horizontal direction was found to be smaller than that in vertical [35]. In addition, the anisotropy of stratospheric turbulence was revealed in [36], with the help of a two-component (isotropic/anisotropic) power spectrum. The validity of such a spectrum was verified by balloon-borne experiments showing that the major contribution to scintillation comes from the anisotropic component.  √ π n 2 L 3 R ox /3R oj , j = y, z, respectively. Here, n 2 is the root-mean square value of the refractive index fluctuations, L is the range, and R oj , j = x, y, z are constants on the order of outer scale. Reprinted with permission from Reference [32].
Recently, several experimental campaigns aimed for extending the results in [32] to the assessment of turbulent air anisotropy in all directions, not solely horizontal and vertical, while doing so in the actual atmosphere. The field measurements of the boundary-layer turbulent anisotropy were carried in Wang et al. [37] on the grounds of the University of Miami, by means of the two-point intensity-intensity correlation function of a nearly spherical wave, i.e., a divergent Gaussian beam (over 200 m path, up to 2 m from the grassy field). The measurements were taken along three horizontal links (differing slightly with respect to effects of the wind penetration) at three heights from the ground. The intensityintensity correlation function (each recorded pixel correlated with the center of the beam) revealed an elliptic form inclined at an angle of about 30 degrees with respect to horizontal. Both the ellipse's eccentricity and orientation angle were shown to be link-and heightdependent. Figure 3 (top row) shows the results of (sequential) measurements carried out for the same path at three heights. The ellipse is the most stretched and rotated for the path closest to the ground. Figure 3 (bottom row) shows the results for three different links taken over the same (middle) height. The ellipse in the two-point intensity-intensity correlations can be directly related to the mean-path averaged anisotropy ellipse of the turbulent eddies.
A similar measurement campaign obtained the intensity covariance function of a plane wave propagating in the vicinity of a runway (2 km path, 2 m above the ground) [38]. As Figure 4 illustrates, a more precise anisotropy structure was revealed in which both positive and negative correlations are seen. The orientation angle was shown to vary up to 90 • to the vertical, depending on the path length and meteorological parameters. The data were well matched with the suggested theoretical model. It was also confirmed in Reference [39] that the ellipse's parameters may have the signature of seasons and time of the day and night. A field measurement of the two-point optical field correlation (mutual coherence) function and the wave structure function (sub 100 m paths, 1.5 m from ground) were also recently carried out at the University of Maryland facilities by means of a plenoptic sensor in Reference [40]. The elliptic shapes, having a substantial orientation angle with respect to the horizontal direction, were also captured in these statistics. While the elliptic shape appearing in various light statistics can be readily explained with the ellipsoid-like eddy correlation function of the turbulent eddies, the qualitative and quantitative justification for the varying orientation angle remains obscure.

Non-Kolmogorov Power Law
In addition to anisotropy effects, the Kolmogorov power law (11/3) of the 3D power spectrum of the refractive index in the inertial range of scales might be violated. In the 1990s, Kyrazis et al. measured non-Kolmogorov turbulence in velocity and temperature statistics in the upper troposphere and lower stratosphere [41,42] and found that substantial qualitative discrepancies are possible between them. Such discrepancies could lead to optical power spectra, with the exponent differing from the Kolmogorov value. The non-Kolmogorov turbulence was also recently measured in an experiment involving urban paths [43].

Presence of Refractive Index Gradients
Yet another effect pertinent to a non-classic turbulent regime is that due to the presence in the same region of space of the refractive index fluctuations and the constant (or very slowly varying) refractive index gradients. The refractive index is typically in the inverse dependence with the average temperature. For example, over a hot ground surface, the average refractive index increases with height. Therefore, an image received by an observer would be shifted down from its actual direction and possibly be overlapped with the direct line-of-sight image. In addition, the turbulent air may add image jitter and diffraction. Recent measurements of this turbulent mirage phenomenon [44] indicate that over the 13-15 km path, the vertical image shift might reach several meters. See also [45,46] for related analytical and numerical calculations.

Deep Turbulence Effects
While performing wave-optics simulation in the strong focusing regime of turbulence, Lachinova and Vorontsov [47] discovered a significant mismatch between numerically estimated and theoretically predicted values of the optical wave scintillation index (normalized variance of fluctuating intensity). The authors showed in their analysis that the reason for this disagreement is due to the irregular appearance of giant irradiance spikes with the amplitudes several times exceeding the diffraction-limited intensity. These spikes would emerge spontaneously due to random formation of focusing lenses extended along the propagation path. Such lenses "trap" irradiance speckles of a suitable size; hence, the spikes can propagate over distances of several kilometers. The statistical analysis of probability of giant spikes appearance provides the distance ranges where these spikes are most likely to be observed. It was also shown that, if compared to the long-exposure (mean) intensity, the spikes whose amplitudes exceed the mean value by a factor of several tens are observed, with the probability being nearly independent from the propagation distance. The probability of spikes' occurrence obtained in the numerical simulations were compared to the theoretical estimations based on the log-normal probability distribution. According to the wave-optics simulations, the giant spikes exceeding the diffraction-limited intensity value by a factor of five or more emerge 10-to 20-fold more frequently than the theory predicts.

Non-Classic Turbulence Emulators
In the laboratory conditions, simulation of a non-classic air turbulence with the controllable anisotropy ratio and the non-Kolmogorov exponent can be achieved by mixing the hot and cold air streams flowing from the opposite directions with different speeds and passing through a fine mesh [48] (see Figure 5). The power law exponent and the anisotropy ratio can then be related to the beam wander analysis of a sufficiently narrow laser beam propagating through mixed air in a direction orthogonal to the air flows. The experimental campaigns indicated that, in fact, achieving the classic turbulence regime precisely is not an easy task, and it does not generally hold for equal temperatures and velocities of the two air streams. Figure 6 shows data for the Hurst exponents H in horizontal and vertical directions in the case of unequal air flow speeds with (top) hot air moving faster and (bottom) cold air moving faster. The Hurst exponent provides an alternative measure of self-similarity (long-term memory) of a random process, calculated directly from the time series. It is limited to interval 0 < H < 1 and characterizes positive/negative correlation for H < 0.5/H > 0.5. The dashed line indicates isotropic turbulence with H = 5/6, which corresponds to the Kolmogorov case. It is evident that in the former case, H is more consistent in both directions and approaches the value of 5/6 for sufficiently high air speeds, while in the later case, the anisotropy is more pronounced but is independent from the air speed. The case of equal air speeds resembles the top subfigure but with more evident tendency to H = 5/6 at high air velocities (not shown).

Theoretical Modeling of Optical Turbulence
In the previous section, the motivation for the thorough analytic modeling of nonclassic optical turbulence was presented. Here, we follow the well-known and recent theories capable of characterizing optical turbulence in various regimes. In particular, we first review the Obukhov-Corrsin theory, adjusting the Kolmogorov spectrum in the neardissipation regime and then outlining the theories for non-Kolmogorov, anisotropic turbulence (in inertial range) followed by jet-stream turbulence and coherent turbulence.

Obukhov-Corrsin Theory
We begin by discussing the Obukhov-Corrsin theory, suitable for characterization of the 3D atmospheric, boundary-free, homogeneous and isotropic turbulence, for which Re can reach values on the order 10 8 -10 9 . Additionally, the Péclet number must be sufficiently high, and the Mach number (both dimensionless) must be subsonic. The Péclet number (for temperature) is defined as a ratio of the rate of advection to the rate of diffusion as follows: and the Mach number is a characteristic of the fluid velocity u relative to the speed of sound u s in the medium as follows: For Ma > 1. 35 and Ma < 0.65 the flow is characterized as supersonic and subsonic, respectively. The fluctuations of the refractive index are conventionally treated as a random process with stationary increments. It is a possibility of involving such a random process since the average value of the refractive index is a slowly-varying function of time. Then, it is sufficient to characterize the refractive-index by its two moments-the average value and the autocovariance function as follows: Here, subscript M stands for the statistical average over the realizations of turbulence, and r is the 3D position vector. If turbulence is homogeneous, using the Wiener-Khinchin theorem, it is also possible to characterize the classic fluctuations in the spatial frequency domain via the power spectrum: r d = r 1 − r 2 and κ κ κ = κ xx + κ yŷ + κ zẑ being the difference vector and the 3D vector of spatial frequencies, respectively. Further, if turbulence is also isotropic, its 3D spectrum can also be written as Φ n (κ), where κ = |κ κ κ|. Based on the Kolmogorov-Obukhov theory for the turbulent velocity field, the power spectrum for any advected scalar field can be determined by the Obukhov-Corrsin theory [49,50], according to which the equilibrium turbulent regime is split into three regimes: inertial-convective, viscous-convective and viscous-diffusive, depending on the participating scales. The inertial-convective regime corresponds to scales between the outer scale L 0 (the largest eddy size) and the following: In cases when the Prandtl number is sufficiently small, the inertial-convective regime can be directly followed by the inertialdiffusive regime. The general form of the 3D power spectrum temperature or for refractive index fluctuations is as follows: where C 0 = 0.72 is the Obukhov-Corrsin constant, χ j is the variance dissipation rate, having units [K 2 s −1 ] for temperature, and [s −1 ] for the refractive index. Function g is a constant in the inertial-convective range function (hence the spectrum is of Kolmogorov −11/3 type), and then it increases to a maximum in the viscous-convective regime and then decreases to zero in the viscous-diffusive regime. Function g depends on Pr. In the inertial-convective regime, the 3D power spectrum is described by Kolmogorov's power law [51] as follows: where C 2 T [m −2/3 ] is the temperature structure parameter as follows: where Γ denotes the Gamma function. For higher spatial frequencies, the spectrum involves function g T : It is a convention in atmospheric propagation studies to use inner scale l 0 instead of η K . The two quantities are proportional: reducing to l 0 = 7.42η K for atmospheric case. The widely used models for function g l (κl 0 ) = g T (κη K ) were suggested in [52] as follows: and in [53] g l (κl 0 ) = exp(1.1090l 0 κ) While Tatarskii's model gives a smooth decay, the Freilich model predicts the occurrence of a bump as spatial frequencies transition from inertial-convective to viscousdiffusing regime (see Figure 7). Other models (e.g., [54]) have been suggested but were later shown to somewhat undershoot/overshoot a constraint that must be imposed on the power spectrum to be consistent with the first principles of thermodynamics [55]. More generally, in the 1978 papers of Hill [56,57], several analytical and numerical models were suggested for more precise modeling of power spectra that can be directly related to the thermodynamic state of the air. For instance, in the maritime atmosphere, the power spectrum may have to be substantially modified at high spatial frequencies, as compared to the over-the-ground turbulence, as was revealed in Reference [58].

Non-Kolmogorov Turbulence
It was shown long time ago by Bolgiano that turbulence does not always follow the Kolmogorov power spectrum model, even in the inertial range [59]. This is due to the fact that the abstraction of turbulent energy by buoyancy forces leads to a sharp decrease in the rate of energy transfer with the wave number. As a consequence, the kinetic energy is transported across the spectrum with a non-uniform rate, while decreasing with increasing the wave number.
Analytic models for the isotropic, non-Kolmogorov power spectrum were used since the 1990s (e.g., [60]). A comprehensive model developed in Reference [61], including the small-and the large-scale cut-offs, has the following form: where C 2 n (z) has units m 3−α and Here, the small and the large spatial frequency cut-offs are as follows: with The major effect of varying α on the turbulence structure is in the weight distribution among the participating turbulent scales: for smaller values of α, more energy is contained in smaller scales, and hence, the turbulence acquires a finer, more granular structure; for larger values of α, more energy is attributed to larger scales, and the turbulence fluctuations takes a cruder structure.
It was suggested in Reference [62] that the non-Kolmogorov power spectrum has a dependence on height, such as z, above the ground, i.e., the following: where and the dependence α = α(z) is obtained by combining the Kolmogorov power spectrum (α 1 = 11/3, A = 0.033) in the boundary layer (up to ≈2 km above the ground), power law (α 2 = 10/3, A = 0.015) for helical turbulence in the troposphere (≈2-8 km) and the power law (α 3 = 5, A = 0.0024) pertinent to the stratosphere (above ≈8 km) as follows: Here, H 1 and H 2 are the adjustable separation altitudes, and b 1 and b 2 are the suitable regime matching coefficients. The corresponding dependence of C 2 n on height z was also discussed in [62] (see also work [63,64] by the same authors). Figure 8 shows the α-height dependence for two sets of parameters b 1 and b 2 .
Some other modeling aspects and representations relating to non-Kolmogorov turbulence were also discussed in [65][66][67].

Anisotropic Turbulence
In the most general case, the turbulent anisotropy can be modeled with the ellipsoidlike correlation function (or the power spectrum) discussed in Section 2.1. As was illustrated above with the optical measurements, such an ellipsoid can be arbitrarily oriented with respect to the ground. However, for analytical modeling, the ellipsoid with the semi-axis describing a vertical and two mutually orthogonal horizontal directions is typically used, in which the vertical semi-axis is shorter than two equal horizontal semi-axes. In this case, a typical turbulent eddy can be viewed as a horizontally oriented "crepe". Such a model was applied to two light propagation cases: (A) vertical path (see Figure 9) and (B) horizontal path (see Figure 10). See also [68] for qualitative theory of anisotropic irregularities.  In case (A), the power spectrum model accounts for the turbulent eddy symmetry axis placed along the optical axis, for example, z [69]. In addition, since the deviation of the power law from the Kolmogorov value of 11/3 is also regarded to be a consequence of stratification and anisotropy [35], these two phenomena can be expressed together as follows: where Γ stands for the Gamma function and C 2 n is the generalized refractive index structure parameter with units m 3−α . The anizotropic factor µ z is related to the ellipse eccentricity in the z (vertical) direction.
In case (B), the power spectrum describes the axis of eddy symmetry being orthogonal to the optical axis z [70,71]: Here, anisotropic factors µ x and µ y define the eccentricity of the ellipse in the x-y planes, i.e., the planes orthogonal to the direction of propagation. In these spectra, the inner scale and the bump at high spatial frequencies are not included but can also be modeled in.
The power spectrum model of Reference [71] assumes the same degree of anisotropy for all turbulent scales in the inertial sub-range, according to Kolmogorov theory [13]. Toselli introduced an extension of this model to the situations when anisotropy may be scale-dependent [72]. In this model, the eddies on the scale comparable to the inner scale are spherically symmetric and gradually increase in the anisotropy in the vertical direction as the outer scale is approached. Other elaborate models for treating anisotropy were also suggested, e.g., in [73]).

Jet-Stream Turbulence
A comprehensive model for turbulence spectrum having anisotropy along all three Cartesian directions for description of a aerojet stream was overviewed in Reference [74], largely based on experimental work by V. S. Sirazetdinov, e.g., in [75]: where measurements return C 2 n on the order of 10 −9 m −2/3 and the longitudinal anisotropy factor Q ≈ 6.

Coherent Turbulence
A number of theoretical models have been proposed for accounting for turbulence containing coherence structures, which may be present in the close proximity of hard boundaries (see References [76,77] and references wherein). Such coherent structures appear because of the temperature or pressure gradients, and are formed as spatial vortices. They have a longer lifespan as compared with classic turbulence eddies. The coherent structures also undergo the direct energy cascade process, being deterministic, unlike the classic turbulent cascade. A typical situation in which a solitary coherent eddy would be generated is an obstacle placed across the air flow. The associated power spectrum has the following form: Here, ν = 1/3 and ν = 5/6 correspond to Kolmogorov and coherent turbulence.

SLM/DMD Benchtop Simulations
Computer simulations of atmospheric optical turbulence and its effects on light beams have been popular since the 1990s (see Reference [78] and references wherein). They allowed to make predictions about the trends in optical beam evolution through a variety of turbulent regimes, both classic and non-classic. An interesting alternative to purely digital simulations are the laboratory benchtop systems, using real light as a source and a spatial light modulator (SLM) or a digital mirror device (DMD) as a simulator of a thin phase screen properly tuned to mimic an extended turbulence path [79]. The non-classic (anisotropic) turbulent thin phase screens were developed and studied in depth, e.g., in [80]. Unlike in computer simulations, having the possibility of using a large set of screens distributed along the propagation path, the SLM/DMD based simulations are limited to one, or possibly two to three screens because of the substantial power loss and the need for synchronization. The fundamental difference between the two devices is their spatial modulation quality and the refresh rate. A quality SLM (with a high fill factor), being a relatively slow device (several hundred Hz), can provide a fine spatial profile, while a DMD, attaining the kHz rates, produces a crude spatial modulation. In addition, at very high refresh rates, the DMD can only operate in a binary mode.
In order to explore the effects of non-Kolmogorov, isotropic and anisotropic turbulence along the vertical paths the laboratory experiment with a liquid crystal, nematic SLM was carried out [81]. The wave optics simulation (WOS) [82] was adopted to the SLM dimensions and the propagation path between the SLM and the camera. It was revealed that for a larger ratio of anisotropic factors, the effect of turbulence on the intensity fluctuations diminishes, i.e., anisotropy acts as a factor to power the spectrum strength. At the same rate, the largest effects on the scintillation index were observed for power law α in the range from 3.1 to 3.2, depending on the applied anisotropy. We remark that the scintillation index shows a peak in the mentioned range if the chosen length unit is the meter [83].
The same numerical procedure and the laboratory arrangement was also used for illustration of the effects of anisotropic turbulence in the horizontal scenario [84]. The anisotropic turbulence was shown to stretch the initially circular beam profile to an elliptical shape oriented along the vertical direction. This is implied by the fact that the weaker/stronger refractive index fluctuations in the horizontal/vertical directions lead to smaller/larger turbulence-induced diffraction: for the fixed ratio of anisotropic factors, µ x and µ y with the largest eccentricity corresponding to α ≈ 3.2 (for fixed ratio of anisotropic factors). Typical phase screens used in these measurement campaigns are given in Figure 11. The single SLM arrangement used in References [81,84] can be also modified to include the second SLM which pre-modulates laser light before sending it through the SLM with the turbulent screens [85]. In this case, a sufficient distance between the two devices must be set to send the pre-modulated light into the far-zone. The iris before the second SLM is used to pass the first-order of the first SLM. The optical setup is controlled from three personal computers: for the two (not synchronized) SLMs and for the camera. Figure 12 illustrates the optical setup and lists examples of an average intensity of the frame-like partially-coherent beams produced by the first SLM, set to interact with non-Kolmogorov turbulence of different power laws, set on the second SLM. Just like for the laser beams, the strongest effect of turbulence is seen for α = 3.2, provided that the chosen length unit is meter.

Theoretical Predictions for Non-Classic Turbulence-Light Interactions
In this section, we will overview the extended Huygens-Fresnel (EHF) method for characterizing the interactions between an optical beam and atmospheric turbulence described with non-Kolmogorov and anisotropic power spectra on the order of the secondorder statistical moments of the field.
Several important second-order statistics can be generally deduced from the knowledge of the cross-spectral density (CSD) matrix [86] of the beam. EHF is currently the most popular method used for predictions of CSD propagation of beams radiated by sources with any spectral composition, degree of coherence and polarization state. For a beam originating in the source plane z = 0 and propagating close to the positive z direction, the EHF integral relates the components of the CSD matrix in the source ← → W (0) and in the field ← → W by the following law: where r α (α = 1, 2) are the 2D source points (z = 0), r α (α = 1, 2) are the 3D field points, ω is the angular frequency, and the integration is taken twice over the source plane. Additionally, i and j refer to the mutually orthogonal components x and y of the electric field, transverse to optical axis z. Propagator K, being the correlation of two spherical wave Green's functions, takes for the 3D power spectrum Φ n (κ κ κ) the following form [87]: Here, r ⊥ d = r ⊥ 2 − r ⊥ 1 , r d = r 2 − r 1 and κ κ κ ⊥ = (κ x , κ y , 0). The first line of propagator K in Equation (40) corresponds to free-space diffraction, while the rest attributes to the effects of turbulence. To arrive at propagation law (39), the Markov approximation was employed, under which the correlation in the forward directions is modeled as the delta-function.
Under the assumption of isotropic turbulence, the turbulence-related part of K reduces further to the following expression: where κ ⊥ = |κ κ κ ⊥ | and J 0 are the zero-order Bessel functions of the first kind. The two-fold integral can be either solved numerically or approximated for small values of the Bessel function argument. The anisotropic spectrum can also be first converted to an isotropic spectrum, and then the procedures above can be applied. For example, in the case of a horizontal link, the transformation can be done by implementing the change of variables µ α κ α → κ α , α = x, y and considering new power spectrum depending only on function |κ| ⊥ in the new variables.
From the knowledge of the field CSD matrix components W ij , the evolution of the spectral density S and the degree of coherence µ can be found as follows: and η(r 1 , r 2 ; ω) = Tr ← → W (r 1 , r 2 ; ω) S(r 1 ; ω) S(r 2 ; ω) .
A variety of polarization properties can also be determined the most important of which is the degree of polarization [86]: Here, Tr and Det stand for matrix trace and determinant.
Second-order statistics of stationary optical beams radiated by a large number of sources with various spectral compositions, spectral density distributions as well as coherence and polarization states were theoretically analyzed in interactions with non-Kolmogorov turbulence (e.g., [88][89][90][91][92][93]). The beams discussed the first three of these references, referring to the broad class of Gaussian Schell-Model (GSM) beams (either vectorial or treated under scalar approximation). These beams are the extensions of the Gaussian laser beams to any coherence state (from coherent to incoherent) and any polarization state (from polarized to unpolarized). They can either be considered at a fixed frequency ω [88,91], or a model spectral composition can be employed [90]. Figure 13 illustrates the average spreading of a GSM beam with growing propagation distance from the source, defined as follows: for several values of α. Figure 14 illustrates the changes in the degree of coherence η and the degree of polarization P as the GSM propagates in non-Kolmogorov turbulence [90]. Interaction of the electromagnetic GSM (EM GSM) beams with anisotropic turbulence along a vertical path (up-link) was examined in Reference [94]. The refractive-index structure constant profile was included in accordance with the well known Hafnagel-Valley model [95]. Figure 15 presents the behavior of the spectral density of a coherent (A) and partially coherent (B) EM GSM beam in the up-link configuration in the presence of the anisotropic, non-Kolmogorov turbulence with α = 3.5 and several values of the effective anisotropy parameter ζ ∼ µ −1 z . For ζ = 1, the turbulence reduces to isotropic, while for ζ > 1, it accounts for anisotropic structure. In the vertical propagation case, the anisotropy only acts as a strength modifier of turbulence. We also note that for a coherent beam, its impact is larger than in the case of a partially coherent beam: in the latter case, most of the diffraction is caused by source correlation at any rate.
The evolution of the EM GSM along a horizontal path in an anisotropic turbulence was considered in Reference [96]. Figure 16 presents the transverse cross-sections of three EM GSM beams with different source coherence states (top row corresponds to more coherent beam) on passing in such turbulence at different ranges. While a more coherent beam acquires an elliptical profile (vertically stretched) on passing at a sufficient distance from the source, the less coherent beam retains its initial circular cross-section, proving to be more resilient to turbulence fluctuations. A similar analysis, performed on the degree of coherence of such beams, leads to a similar set of cross-sections but the ellipses are stretched along the horizontal direction [96].  Light propagation analysis through turbulence with the jet-stream spectrum slightly modified for integral convergence was performed in [97]. In addition, the behavior of higher-order statistics of light in non-classic turbulence was also studied by other methods, e.g., the Rytov method. We will discuss some of these results below on considering the applications. The evolution of the second-order statistics of non-stationary (pulsed) fields in non-Kolmogorov turbulence was treated in [98,99].

Wireless Optical Communications
Wireless optical communication (WOC) systems using laser light are attractive because of high data rates (due to extremely wide bandwidth) and enhanced security (supported by high laser directionality). However, WOC systems suffer from atmospheric turbulenceinduced signal degradation. Non-classic turbulent regimes can also affect WOC in a number of ways. Toselli and Korotkova used the extended non-Kolmogorov anisotropic power spectrum introduced in [71] to investigate the performance of a WOC system (with a Gaussian beam used as the information carrier) embedded into anisotropic turbulence [100]. Specifically, the scintillation index/flux with aperture averaging, the probability of fade, the mean signal-to-noise ratio (SNR) and the bit-error rate (BER) for the on-off key modulation were analyzed for a specific case of horizontal propagation path at high altitude.

Aperture-Averaged Scintillation of a Gaussian Beam
The intensity fluctuations measured at the focal plane of the receiving telescope are conventionally characterized by the aperture-averaged scintillation flux [100], generally defined as follows: with I being the instantaneous optical intensity (see [71] containing the well-known expressions for σ I derived for basic types of waves). This parameter is crucial since it is required for calculation of the probability of fade, the SNR and the BER in the direct detection systems. The complete theoretical analysis of the scintillation flux in non-classic turbulence is given in [100] (see, in particular, Equation (14)). Here, we display in Figure 17 the scintillation flux with the fixed structure function at the Fresnel distance, for six pairs of anisotropic factors along the x and y axes (shown as z x eff and z y eff in the figure legend). It is evident that at higher values of anisotropy factors, the scintillation flux is reduced at any power law in the range from 3 to 4.

Probability of Fade
For a given probability density function (PDF) of the irradiance fluctuations, the probability of fade describes the percentage of time, and the irradiance of the received signal is below some prescribed threshold value, for example, I t . Hence, the probability of fade as a function of a set threshold level is defined by the cumulative probability as follows: where P I (I) is the PDF of the fluctuating instantaneous irradiance. In weak turbulence, the log-normal PDF model is most frequently used for description of the optical intensity fluctuations, leading to the following expression for the probability of fade [95]: where Erf(x) is the error function. Here, threshold (fade) parameter F t represents the intensity margin in decibels (dB) from the threshold I t , which is usually the sensitivity of the receiver. On using the scintillation flux (from [100]) in Equation (48), one can deduce the probability of fade as a function of slope α for a fixed fade threshold parameter (say F t = 6 dB) and fixed structure function at the Fresnel distance, in a particular horizontal case scenario and for different values of anisotropy parameters along the x and y axes (shown as z x eff and z y eff in the legend). The curves in Figure 18 indicate that the decrease in scintillation due to anisotropy reduces the probability of fade.

Mean-Signal-to-Noise Ratio
Due to turbulent fluctuations, the received irradiance must be treated as a random variable, over long detection intervals. Hence, in the case of a shot-noise limited system and under the assumption of sufficiently small beam spreading due to turbulence, the mean SNR at the output of the detector SNR assumes the following general form [95]: where the scintillation flux is developed in [100] and SNR 0 is the free-space SNR. We show in Figure 19 the mean SNR, SNR , keeping the structure function at the Fresnel distance fixed, as shown in [100] and previously proposed in [83]. The two anisoropic factors are kept fixed (shown as z x eff and z y eff in the figure title). The impact of anisotropy on the mean SNR is noticeable: in general, the decrease in scintillation due to anisotropy leads to an increase in the SNR, implying better performance.

Mean Bit-Error Rate
The probability of error of a beam wave after propagation in the atmospheric turbulence is the conditional probability that must be averaged over the PDF of the random received signal to determine the unconditional mean BER. In terms of a normalized signal with unit mean, if a direct detection scheme is used, this leads to the following expression [95]: where Erfc is the complimentary error function and the PDF of the instantaneous intensity is assumed to be log-normal distributed with the unit mean. In Figure 20, we show the mean BER for different values of power law α as a function of the mean SNR (in dB), for several sets of anisotropy factors (shown as z x eff and z y eff in the legend) and a given value of the structure function (kept fixed) [100]. The effect of anisotropy on the mean bit-error rate is well visible: as for the probability of fade and the signal-to-noise ratio, the decrease in scintillation, due to anisotropy leads to a decrease in the bit-error rate. Toselli and Gladysz also showed that the removal by adaptive optics of several Zernike modes from the turbulence-affected phase can be very effective in reducing scintillation in non-classic turbulence [101].

LIDARs
Light detection and ranging systems (LIDARs) are currently used in various applications from meteorology to headlight sensing. In these systems, light must travel through the same medium twice, either through exactly the same path (mono-static case) or along slightly different paths (bi-static case). In the former case, the phase conjugation phenomenon considerably complicates the analysis, producing intensity redistribution in the region around the optical axis. This phenomenon is known as enhanced back-scatter (EBS), and is positive for retro-reflector targets and negative for flat mirror targets. While double-passage propagation of laser beams were thoroughly investigated for Kolmogorov turbulence [95,102] a more complex case of non-Kolmogorov turbulence was not fully explored until recently. The seminal analysis of LIDARs in weak non-Kolmogorov turbulence was made in [103], revealing how a power law different from Kolmogorov's impacts the average intensity distribution, the long-term spread and the scintillation index of an optical beam after it passes through turbulence to a (smooth) target surface, interacting with such a target and passing back through the same turbulent channel. Such an analysis was conducted for both mono-static and bi-static configurations. In the case of the mono-static channel, the enhanced backscattering effects were also investigated. Both the analytical results and those based on wave-optics simulation were obtained and compared.
The double-passage problem of a laser beam propagating in the presence of non-Kolmogorov atmospheric turbulence at any turbulence strength was analyzed in Reference [104]. In particular, using the extended Rytov theory, for the horizontal path, the authors theoretically investigated the scintillation index of a Gaussian beam reflected from a small unresolved target in deep turbulence conditions. The authors theoretically showed that different power laws substantially affect the scintillation results; however, they did not observe the scintillation peaks predicted by the theory near the focusing regime, using wave optics simulations. In addition, authors showed that the occurrence of giant intensity spikes previously discussed in Section 2.4 is currently not captured by the theory. The numerical outcomes for the scintillation index in mono-static and bi-static cases, at any turbulence strength, are shown in Figures 21 and 22, respectively. Finally, we show in Figure 23 the EBS factor appearing in the mono-static LIDAR scenario, as a function of the square root of the Rytov variance for several power law values [104].

Imaging Systems
The analysis of imaging systems using random illumination and/or operating through random environments relies on the modulation transfer function (MTF) describing a filter acting on the continuum of spatial frequencies characterizing the system [95]. The MTF can be directly derived from the CSD function, discussed above. In the clear-air (particle-free) atmospheric turbulence, the MTF is defined by the following convolution: where the first term prescribes the properties of the system in the absence of turbulence, and the second term characterizes spatial filtering due to turbulence, its form depending on long or short exposure options. Parameter Ω represents the spatial frequency scaled by wavelength and the collecting aperture diameter. Imaging through various regimes of nonclassic turbulence was discussed in References [60,[105][106][107][108]. The influence of both non-Kolmogorov power spectra and anisotropy on the MTF profiles, and thus, on the quality of the formed images, is illustrated.

Concluding Remarks
We have overviewed the classic and a relatively recent body of work highlighting various manifestations of non-classic optical turbulence in the atmosphere and discussed the ways of characterizing and predicting its effects on propagating light. The main accent was made on phenomena and applications that were investigated by the authors over the last decade, through analytical modeling and computer simulations as well as in lab and field measurement campaigns.
Optical wavefronts are extremely susceptible to turbulent air fluctuations and, hence, present excellent means for the sensing of classic and non-classic turbulence. We have brought together a large number of investigations reported in the literature, which illustrate that the exponent variations of the power spectrum, anisotropy, constant temperature gradients, and other manifestations of non-classic turbulence are imprinted into light statistics and can be used for assessing turbulence's structure and statistics. On the other hand, non-classic turbulent regimes present obstacles for tuning various optical systems operating through the air. As we showed, careful theoretical analysis and simulations can provide an idea of the limits that the non-classic turbulence can set for high-quality operation of WOC, imaging and LIDAR systems. To conclude, we must mention that other natural turbulent media, such as ocean water and soft biological tissues, often exhibit non-classic regimes. Therefore, the summary suggested here may also become useful for further investigations relating to these media.