A New Method of Generating Atmospheric Turbulence with a Liquid Crystal Spatial Light Modulator

1. Introduction

Light traveling from a star, or any point source, will propagate spherically outward. After a long distance, the wavefront, or surface of equal phase, will be flat; as is illustrated in Fig. 1.

Figure 1.

Flat wavefront after a long propagation distance from a point source

When the light begins to propagate through Earth’s atmosphere, the varying index of refraction will alter the optical path, as shown in Fig. 2. The Earth’s atmosphere can be described as a locally homogeneous medium in which its properties vary with respect to temperature, pressure, wind velocities, humidity and many other factors. Also, the Earth’s atmosphere temporally changes in a quasi-random fashion. All of these processes are usually simply refered to as ``atmospheric turbulence’’. The Kolmogorov model of energy distribution in a turbulent medium is a useful statistical model to describe the fluctuation in refractive index due to mostly the humidity and pressure changes. This model was first proposed by a Russian mathematician named Andreï Kolmogorov in 1941 and describes how in a fully turbulent media the kinetic energy of large scale motions is transfered to smaller and smaller scale motions (Kolmogorov, 1941). It is supported by a variety of experimental measurements and is quite widely used in simulations for the propagation of electromagnetic waves through a random medium. The first author to fully describe such phenomena was Tatarski in his textbook ``Wave propagation in a turbulent medium’’ (Tatarski, 1961). The complex and random nature of the Earth’s atmospheric turbulence effect on wave propagation is currently a subject of active research and experimental measurements. Many of the parameters of Earth’s atmospheric turbulence can be, at best, described statistically.

Figure 2.

Propagation of light from a distant source that then passes through the atmosphere

These statistical parameters represent the strength and changeability of the atmospheric turbulence, these conditions are customarily refered to as the ``astronomical seeing’’, as they are widely used for astronomical applications. It is with this statistical information about a certain astronomical site and the specifications of the telescope that an Adaptive Optics (AO) system can be designed to correct the wavefront distortions caused by the atmosphere at that site. As telescopes continue to be manufactured larger and larger, the need for AO is increasing because of the limiting factors caused by atmospheric turbulence. In order to adequately characterize the performance of a particular AO system, an accurate spatial and temporal model of the Earth’s atmosphere is required.

AO is the term used for a class of techniques dealing with the correction of wavefront distortions in an optical system in real time. Some wavefront distortions may include those caused by the atmosphere. Astronomical applications of AO particularly include the correction of atmospheric turbulence for a telescope system. Other possible applications include Free Space Laser Communications, High Energy Laser Applications, and Phase-Correction for Deployable Space-Based Telescopes and Imaging systems. However, prior to deployment, an AO system requires calibration and full characterization in a laboratory environment.

Many techniques are currently being used with AO systems for simulating atmospheric turbulence. Some static components use glass phase screens with holograms etched into them. In addition, it is also important to simulate the temporal transitions of atmospheric turbulence. Some of these methods include the use of a static aberrator, such as a clear piece of plastic or glass etched phase screen, and rotating it. Rotating filter wheels with etched holographic phase screens can simulate temporal transitions, as well. Also, simply using a hot-plate directly under the beam path in an optical system can simulate temporally the atmospheric turbulence.

However, etching holographic phase screens into glass can be quite costly and not very flexible to simulate different atmospheric characteristics. Thus, one would need more than one phase screen. A testbed that simulates atmospheric aberrations far more inexpensively and with greater fidelity and flexibility can be achieved using a Liquid Crystal (LC) Spatial Light Modulator (SLM). This system allows the simulation of atmospheric seeing conditions ranging from very poor to very good and different algorithms may be easily employed on the device for comparison. These simulations can be dynamically generated and modified very quickly and easily.

2. Background

2.1. Brief history of the study of atmospheric turbulence

Ever since Galileo took a first look at the moons of Jupiter through one of the first telescopes, astronomers have strived to understand our universe. Within the last century, telescopes have enabled us to learn about the far reaches of our universe, even the acceleration of the expansion of the universe, itself. The field of building telescopes has been advancing much in recent years. The twin Keck Telescopes on the summit of Hawaii’s dormant Mauna Kea volcano measure 10 meters and are currently the largest optical telescopes in the world. Plans and designs for building 30 and 100 meter optical telescopes are underway. As these telescope apertures continue to grow in diameter, the Earth’s atmosphere degrades the images we try to capture more and more. As Issac Newton said is his book, Optiks in 1717, “… the air through which we look upon the stars is in perpetual tremor; as may be seen by the tremulous motion of shadows cast from high towers, and by the twinkling of the fixed stars…. The only remedy is a most serene and quiet air, such as may perhaps be found on the tops of high mountains above grosser clouds.” It was at this time when we first realized that the Earth’s atmosphere was the major contributor to image quality for ground-based telescopes. The light arriving from a distant object, such as a star, is corrupted by turbulence-induced spatial and temporal fluctuations in the index of refraction of the air.

In 1941, Kolmogorov published his treatise on the statistics of the energy transfer in a turbulent flow of a fluid medium. Tatarskii used this model to develop the theory of electro-magnetic wave propagation through such a turbulent medium. Then, Fried used Tatarskii’s model to introduce measurable parameters that can be used to characterize the strength of the atmospheric turbulence.

The theory of linear systems allows us to understand how a system transforms an input just by defining the characteristic functions of the system itself. Such a characteristic function is represented by a linear operator operating on an impulse function. The characteristic system function is generally called the ``impulse response function’’. Very often, such an operator is the so-called Fourier transform. An imaging system can be approximated by a linear, shift-invariant system over a wide range of applications. The next few sections will explain the use of a Fourier transform in such an optical imaging system and its applications with optical aberrations.

2.2. Brief overview of fourier optics and mathematical definitions

A fantastic tool for the mathematical analysis of many types of phenomena is the Fourier transform. The 2-dimensional Fourier transform of the function g(x,y) is defined as,

$$G\left(f_x,f_y\right)=\mathbb{F}\left\{g\left(x,y\right)\right\}={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}g\left(x,y\right)e^{-j2\pi \left(f_{x}x+f_{y}y\right)}dxdy}}$$ E1

where, for an imaging system, the x-y plane is the entrance pupil and the f _x -f _y plane is the imaging plane. A common representation of the Fourier transform of a function is by the use of lower case for the space domain and upper case for the Fourier transform, or frequency domain. Similarly, the inverse Fourier transform of the function G(f _x ,f _y ) is defined as,

$$g\left(x,y\right)={\mathbb{F}}^{-1}\left\{G\left(f_x,f_y\right)\right\}={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}G\left(f_x,f_y\right)e^{j2\pi \left(f_{x}x+f_{y}y\right)}d{f}_{x}d{f}_y}}$$ E2

There exist various properties of the Fourier transform. The linearity property states that the Fourier transform of the sum of two or more functions is the sum of their individual Fourier transforms and is shown by,

$$\mathbb{F}\left\{ag\left(x,y\right)+bf\left(x,y\right)\right\}=a\mathbb{F}\left\{g\left(x,y\right)\right\}+b\mathbb{F}\left\{f\left(x,y\right)\right\}$$ E3

where a and b are constants. The scaling property states that stretching or skewing of a function in the x-y domain results in skewing or stretching of the Fourier transform, respectively, and is shown by,

$$\mathbb{F}\left\{g\left(ax,by\right)\right\}=\frac{1}{\left|ab\right|}G\left(\frac{f_x}{a},\frac{f_y}{b}\right)$$ E4

where a and b are constants. The shifting property states that the translation of a function in the space domain introduces a linear phase shift in the frequency domain and is shown by,

$$\mathbb{F}\left\{g\left(x-a,y-b\right)\right\}=G\left(f_x,f_y\right)e^{-j2\pi \left(a{f}_x+b{f}_y\right)}$$ E5

where a and b are constants. This property is of particular interest in the mathematical analysis of tip and tilt in an optical system, as it describes horizontal or vertical position in the imaging plane. Parsaval’s Theorem is generally known as a statement for the conservation of energy and is shown as,

$${\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\left|g\left(x,y\right)\right|}^{2}dxdy=}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\left|G\left(f_x,f_y\right)\right|}^{2}d{f}_{x}d{f}_y}}$$ E6

The convolution property states that the convolution of two functions in the space domain is exactly equivalent to the multiplication of the two functions’ Fourier transforms, which is usually a much simpler operation. The convolution of two functions is defined as,

$$g\left(x,y\right)\ast f\left(x,y\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}g\left(\xi ,\eta \right)f\left(\xi -x,\eta -y\right)d\xi d\eta }}$$ E7

The convolution property is shown as,

$$\mathbb{F}\left\{g\left(x,y\right)\ast f\left(x,y\right)\right\}=G\left(f_x,f_y\right)F\left(f_x,f_y\right)$$ E8

A special case of the convolution property is known as the autocorrelation property and is shown as,

$$\mathbb{F}\left\{g\left(x,y\right)\ast g^{*}\left(x,y\right)\right\}={\left|G\left(f_x,f_y\right)\right|}^2$$ E9

where the superscript ^* denotes the complex conjugate of the function g(x,y). The autocorrelation property gives the Power Spectral Density (PSD) of a function and is a useful way to interpret a spatial function’s frequency content. The square of the magnitude of the G(f _x ,f _y ) function is also referred to as the Point Spread Function (PSF). The PSF is the imaging equivalent of the impulse response function. It is easy to see that the PSF represents the spreading of energy on the output plane of a point source at infinity.

The spatial variation as a function of spatial frequency is described by the Optical Transfer Function (OTF). The OTF is defined as the Fourier transform of the PSF written as,

$$\text{OTF}=\mathbb{F}\left\{{\left|G\left(f_x,f_y\right)\right|}^2\right\}=\mathbb{F}\left\{\text{PSF}\right\}$$ E10

The Modulation Transfer Function (MTF) is the magnitude of the OTF and is written as,

$$\text{MTF}=\left|\mathbb{F}\left\{{\left|G\left(f_x,f_y\right)\right|}^2\right\}\right|=\left|\mathbb{F}\left\{\text{PSF}\right\}\right|$$ E11

Two common aperture geometries, or pupil functions, that will be discussed are the rectangular and circular apertures. The rectangular aperture is defined as,

$$rect\left(\frac{x}{k},\frac{y}{l}\right)=\{\begin{array}{c}1\\ 0\end{array}\begin{array}{c}\left|x\right|\le \raisebox{1ex}{$k$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\text{ and }\left|y\right|\le \raisebox{1ex}{$l$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\\ \text{ otherwise }\end{array}$$ E12

where k and l are positive constants that refer to the length and width of the aperture, respectively. The circular aperture is defined as,

$$circ\left(\frac{\rho }{l}\right)=\{\begin{array}{c}1\\ 0\end{array}\begin{array}{c}\rho \le l\text{ and }\rho =\sqrt{x^2+y^2}\\ \text{ otherwise }\end{array}$$ E13

where l is a positive constant referring to the radius of the aperture. These pupil functions become of great use when analyzing an imaging system with these apertures. For the purposes of this discussion, a circular aperture will be considered as it is of particular use with Zernike polynomials and Karhunen-Loeve polynomials which will be discussed later.

In order to include the effects of aberrations, it is useful to introduce the concept of a “generalized pupil function”. Such a function is complex in nature and the argument of the imaginary exponential is a function that represents the optical phase aberrations by,

$$\mathbb{P}\left(x,y\right)=P\left(x,y\right)e^{j\frac{2\pi }{\lambda }W\left(x,y\right)}$$ E14

where P(x,y) = circ(ρ), λ is the wavelength, and W(x,y) is the effective path length error, or error in the wavefront. It is in this wavefront error that atmospheric turbulence induces and degrades image quality of an optical system and induces aberrations. This wavefront error can be induced in an optical system through the use of a LC SLM. The next several sections will describe methods of simulating atmospheric turbulence in an optical system and introduce the new method of simulating atmospheric turbulence developed at the Naval Research Laboratory.

2.3. Optical aberrations as Zernike polynomials

The primary goal of AO is to correct an aberrated, or distorted, wavefront. A wavefront with aberrations can be described by the sum of an orthonormal set of polynomials, of which there are many. One specific set is the so called Zernike polynomials, Z _i (ρ,θ), and they are given by,

$$Z_i\left(\rho ,\theta \right)=\{\begin{array}{l}\sqrt{n+1}{R}_n^m\left(\rho \right)\cos \left(\theta \right)\hfill \\ \sqrt{n+1}{R}_n^m\left(\rho \right)\sin \left(\theta \right)\hfill \\ R_n^0\left(\rho \right)\hfill \end{array}\begin{array}{c}\\ \\ \end{array}\begin{array}{c}\text{ for }m\ne 0\text{ and }i\text{ is even}\\ \text{ for }m\ne 0\text{ and }i\text{ is even}\\ \text{for }m=0\end{array}$$ E15

where

$$R_n^m\left(\rho \right)={\displaystyle \sum _{s=0}^{\frac{n-m}{2}}\frac{{\left(-1\right)}^s\left(n-s\right)!}{s!\left(\frac{n+m}{2}-s\right)!\left(\frac{n-m}{2}-s\right)!}{\rho }^{n-2s}}$$ E16

The azimuthal and radial orders of the Zernike polynimials, m and n, respectively, satisfy the conditions that m ≤ n and n-m = even, and i is the Zernike order number (Roggemann & Welsh, 1996). The Zernike polynomials are used because, among other reasons, the first few terms resemble the classical aberrations well known to lens makers. The Zernike order number is related to the azimuthal and radial orders via the numerical pattern in Table 1.

i	n	m	i	n	m	i	n	m	i	n	m
1	0	0	8	3	-1	15	4	-4	22	6	0
2	1	1	9	3	3	16	5	1	23	6	2
3	1	-1	10	3	-3	17	5	-1	24	6	-2
4	2	0	11	4	0	18	5	3	25	6	4
5	2	2	12	4	2	19	5	-3	26	6	-4
6	2	-2	13	4	-2	20	5	5	27	6	6
7	3	1	14	4	4	21	5	-5	28	6	-6

Table 1.

Relationship between Zernike order and azimuthal and radial orders

Zernike polynomials represent aberrations from low to high order with the order number. A wavefront can generally be represented by,

$$\text{Wavefront}\left(\rho ,\theta \right)={\displaystyle \sum _{i=1}^M{a}_i{Z}_i\left(\rho ,\theta \right)}$$ E17

where the a _i ’s are the amplitudes of the aberrations and M is the total number of Zernike orders the wavefront is represented by. This wavefront can be substituted into Equation (14) as it represents the phase in an imaging system.

2.4. Kolmogorov’s statistical model of atmospheric turbulence

The Sun’s heating of land and water masses heat the surrounding air. The buoyancy of air is a function of temperature. So, as the air is heated it expands and begins to rise. As this air rises, the flow becomes turbulent. The index of refraction of air is very sensitive to temperature. Kolmogorov’s model provides a great mathematical foundation for the spatial fluctuations of the index of refraction of the atmosphere. The index of refraction of air is given by,

$$n\left(\overrightarrow{r},t\right)=n_0+n_1\left(\overrightarrow{r},t\right)$$ E18

where $$\overrightarrow{r}$$ is the 3-dimensional space vector, t is time, n ₀ is the average index of refraction, and $$n_1\left(\overrightarrow{r},t\right)$$ is the spatial variation of the index of refraction. For air, we may say n ₀ = 1. At optical wavelengths the dependence of the index of refraction of air upon pressure and temperature is n ₁ = n – 1 = 77.6x10^-6 P/T, where P is in millibars and T is in Kelvin. The index of refraction for air can now be given as,

$$n\left(P,T\right)=1+\frac{77.6\times {10}^{-6}P}{T}$$ E19

Differentiating the index of refraction with respect to temperature gives,

$$\frac{\partial }{\partial T}n\left(P,T\right)=-\frac{77.6\times {10}^{-6}P}{T^2}$$ E20

From Equation (20), we can see that the change in index of refraction with respect to temperature cannot be ignored (Roggemann & Welsh, 1996). These slight variances of temperature, of which the atmosphere constantly has many, will affect the index of refraction enough to affect the resolution of an imaging system.

As light begins to propagate through Earth’s atmosphere, the varying index of refraction will alter the optical path slightly. To a fairly good approximation, the temperture and pressure can be treated as random variables. Unfortunately, because of the apparent random nature of Earth’s atmosphere, it can at best be described statistically. It is with this statistical information about a certain astronomical site and the specifications of the telescope that an adaptive optics system can be designed to correct the wavefront distortions caused by the atmosphere at that site.

The quantity $$C_n^2$$ is called the structure constant of the index of refraction fluctuations with units of m^−2/3 (Roggemann & Welsh, 1996), it is a measurable quantity that indicates the strength of turbulence with altitude in the atmosphere. The value $$C_n^2$$ can vary from ~10⁻¹⁷ m^−2/3 or less and ~10⁻¹³ m^−2/3 or more in weak and strong conditions, respectively (Andrews, 2004). $$C_n^2$$ can have peak values during midday, have near constant values at night and minimum values near sunrise and sunset. These minimum values’ occurrence at sunrise and sunset is known as the diurnal cycle.

Figure 3.

Plots of the Hufnagel-Valley, SLC-Day, SLC-Night, and Greenwood models for $$C_n^2$$ with respect to altitude.

Some commonly accepted models of $$C_n^2\left(h\right)$$ as functions height are the Hufnagel-Valley, SLC-Day, SLC-Night and Greenwood models. The Hufnagel-Valley model is written as,

$$C_n^2\left(h\right)=0.00594{\left(\raisebox{1ex}{$v_w$}\!\left/ \!\raisebox{-1ex}{$27$}\right.\right)}^2{\left({10}^{-5}h\right)}^{10}{e}^{\raisebox{1ex}{$-h$}\!\left/ \!\raisebox{-1ex}{$1000$}\right.}+2.7\times {10}^{-16}{e}^{\raisebox{1ex}{$-h$}\!\left/ \!\raisebox{-1ex}{$1500$}\right.}+C_n^2\left(0\right)e^{\raisebox{1ex}{$-h$}\!\left/ \!\raisebox{-1ex}{$100$}\right.}$$ E21

where v _w is the rms wind speed and $$C_n^2\left(0\right)$$ is the ground-level value of the structure constant of the index of refraction. The SLC-Day model is written as,

$$C_n^2\left(h\right)=\{\begin{array}{l}1.7\times {10}^{-14}\hfill \\ 3.13\times {10}^{-13}{h}^{-1.05}\hfill \\ 1.3\times {10}^{-15}\hfill \\ 8.87\times {10}^{-7}{h}^{-3}\hfill \\ 2.0\times {10}^{-16}{h}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\hfill \end{array}\begin{array}{c}\\ \\ \\ \\ \end{array}\begin{array}{c}\\ \\ \\ \\ \end{array}\begin{array}{c}\\ \\ \\ \\ \end{array}\begin{array}{c}\\ \\ \\ \\ \end{array}\begin{array}{c}0h18.5\\ 18.5h240\\ 240h880\\ 880h7200\\ 7200h20000\end{array}$$ E22

The SLC-Night model is written as,

$$C_n^2\left(h\right)=\{\begin{array}{l}8.4\times {10}^{-15}\hfill \\ 2.87\times {10}^{-12}{h}^{-2}\hfill \\ 2.5\times {10}^{-16}\hfill \\ 8.87\times {10}^{-7}{h}^{-3}\hfill \\ 2.0\times {10}^{-16}{h}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\hfill \end{array}\begin{array}{c}\\ \\ \\ \\ \end{array}\begin{array}{c}\\ \\ \\ \\ \end{array}\begin{array}{c}\\ \\ \\ \\ \end{array}\begin{array}{c}\\ \\ \\ \\ \end{array}\begin{array}{c}0h18.5\\ 18.5h110\\ 110h1500\\ 1500h7200\\ 7200h20000\end{array}$$ E23

The Greenwood model is written as,

$$C_n^2\left(h\right)=\left[2.2\times {10}^{-13}{\left(h+10\right)}^{-1.3}+4.3\times {10}^{-17}\right]e^{-\raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$4000$}\right.}$$ E24

In each of these models, h may be replaced by $$\frac{h}{\cos \left({\theta }_z\right)}$$ if the optical path is not vertical, or at zenith, and θ _z is the angle away from zenith.

2.5. Fried and Noll’s model of turbulence

The fact that a wavefront can be expressed as a sum of Zernike polynomials is the basis for Noll’s analysis on how to express the phase distortions due to the atmosphere in terms of Zernike polynomials.

Fried’s parameter, also known as the coherence length of the atmosphere and represented by r ₀ , is a statistical description of the level of atmospheric turbulence at a particular site. Fried’s parameter is given by,

$$r_0={\left[0.423k^2\sec \zeta {\displaystyle \underset{\text{Path}}{\int }{C}_n^2\left(z\right)dz}\right]}^{-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$5$}\right.}$$ E25

where $$k=\frac{2\pi }{\lambda }$$ and λ is the wavelength, ζ is the zenith angle, the Path is from the light source to the telescope’s aperture along the z axis and it is expressed in centimeters. The value of r ₀ ranges from under 5 cm with poor seeing conditions to more than 25 cm with excellent seeing conditions in the visible light spectrum. The coherence length limits a telescope’s resolution such that a large aperture telescope without AO does not provide any better resolution than a telescope with a diameter of r ₀ (Andrews, 2004). In conjunction with r ₀ , another parameter that is important is the isoplanatic angle, θ ₀ , given and approximated by,

$${\theta }_0={\left[2.91k^2{\sec }^{\raisebox{1ex}{$8$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\zeta {\displaystyle \underset{\text{Path}}{\int }{C}_n^2\left(z\right)z^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}dz}\right]}^{-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$5$}\right.}\approx 0.4125r_0$$ E26

and is expressed in milli-arcseconds. The isoplanatic angle describes the maximum angular difference between the paths of two objects in which they should traverse via the same atmosphere. This is illustrated in Fig. 4.

It is also important to remember that the atmosphere is a statistically described random medium that has temporal dependence as well as spatial dependence. One common simplification is to assume that the wind causes the majority of the distortions, temporally. The length of time in which the atmosphere will remain roughly static is represented by τ ₀ and is approximated by,

Figure 4.

Illustration of isoplanatic angle

Zernike Mode	Zernike-Kolmogorov residual error
Tip	Δ1 = 1.0299 (D/r0)5/3
Tilt	Δ2 = 0.5820 (D/r0)5/3
Focus	Δ3 = 0.1340 (D/r0)5/3
Astigmatism X	Δ4 = 0.0111 (D/r0)5/3
Astigmatism Y	Δ5 = 0.0880 (D/r0)5/3
Coma X	Δ6 = 0.0648 (D/r0)5/3
Coma Y	Δ7 = 0.0587 (D/r0)5/3
Trefoil X	Δ8 = 0.0525 (D/r0)5/3
Trefoil Y	Δ9 = 0.0463 (D/r0)5/3
Spherical	Δ10 = 0.0401 (D/r0)5/3
Secondary Astigmatism X	Δ11 = 0.0377 (D/r0)5/3
Secondary Astigmatism Y	Δ12 = 0.0352 (D/r0)5/3
Higher orders (J "/ 12)	ΔJ = 0.2944(D/r0)5/3

Table 2.

Zernike-Kolmogorov residual errors, Δ _J , and their relation to D/r ₀

$$J^{\sqrt{3/2}}$$ E27

where v _w is the average wind speed at ground level, and D is the telescope aperture. The three parameters r ₀ , θ ₀ , and τ ₀ are required to know the limitations and capabilities of a particular site in terms of being able to image objects through the atmosphere.

To make a realization of a wavefront after being distorted by the Earth’s atmosphere, Fried derived Zernike-Kolmogorov residual errors (Fried, 1965, Noll, 1976, Hardy, 1998). The a _i ’s in Equation (17) are calculated from the Zernike-Kolmogorov residual errors, Δ _J , measured through many experimental procedures and calcutated by Fried (Fried, 1965) and by Noll (Noll, 1976) and are given in Table 2. Thus, a realization of atmospheric turbulence can be simulated for different severities of turbulence and for different apertures.

2.6. Frozen Seeing model of atmospheric turbulence

Time dependence of atmospheric turbulence is very complex to simulate and even harder to generate in a laboratory environment. One common and widely-accepted method of simulating temporal effects of atmospheric turbulence is by the use of Frozen Seeing, also known as the Taylor approximation (Roggemann & Welsh, 1996). This approximation assumes that given a realization of a large portion of atmosphere, it drifts across the aperture of interest with a constant velocity determined by local wind conditions, but without any other change, whatsoever (Roddier, 1999). This technique has proved to be a good approximation given the limited capabilities of simulating accurate turbulence conditions in a laboratory environment. For example, a large holographic phase screen can be generated and may be simply moved across an aperture and measurements can then be made. A sample realization of atmospheric turbulence with a ratio of D/r ₀ = 2.25 can be seen in Fig. 5.

Figure 5.

A sample phase screen generated via the Frozen Seeing method

3. New method of generating atmospheric turbulence with temporal dependence

In this next section, a new method of generating atmospheric turbulence is introduced. This method takes into account the temporal and spatial effects of simulating atmospheric turbulence with the thought in mind of being able to use this method in a laboratory with a LC SLM. Some advantages of this method include far less computational constraints than using the Frozen Seeing model in software. In addition, the use of Karhunen-Loeve polynomials is introduced rather than using Zernike polynomials, as they are a statistically independent set of orthonormal polynomials.

3.1. Karhunen-Loeve polynomials

Karhunen-Loeve polynomials are each a sum of Zernike polynomials, however, they have statistically independent coefficients (Roddier, 1999). This is important due to the nature of atmospheric turbulence as described by the Kolmogorov model following Kolmogorov statistics. The Karhunen-Loeve polynomials are given by,

$${\tau }_0={\left[2.91k^2\sec \zeta {\displaystyle \underset{\text{Path}}{\int }{C}_n^2\left(z\right)v_w{}^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}dz}\right]}^{-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$5$}\right.}\approx \frac{0.314r_0}{v_w}{\left(\frac{D}{r_0}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}$$ E28

where the b _p,j matrix is calculated and given by Wang and Markey (Wang & Markey, 1978), and N is the number of Zernike orders the Karhunen-Loeve order j is represented by. Thus, to represent a wavefront, Equation (17) can be rewritten as,

$$K_p\left(\rho ,\theta \right)={\displaystyle \sum _{j=1}^N{b}_{p,j}{Z}_j\left(\rho ,\theta \right)}$$ E29

and now the wavefront is now represented as a sum of Karhunen-Loeve polynomials with the Zernike-Kolmogorov resitual error weights in the a _i ’s.

3.2. Spline technique

Tatarski’s model describes the phase variances to have a Gaussian random distribution (Tatarski, 1961). So, by taking Equation (29) and modifying it such that there is Gaussian random noise factored in gives,

$$\text{Wavefront}\left(\rho ,\theta \right)={\displaystyle \sum _{i=1}^M{a}_i{K}_i\left(\rho ,\theta \right)}$$ E30

where X _i is the amount of noise for the i ^th mode based on a zero-mean unitary Gaussian random distribution and the a _i ’s are the amplitudes of the aberrations calculated from Zernike-Kolmogorov residual errors in Table 2.

The X _i ’s in Equation (30) can be generated by just using randomly generated numbers. But generating a continuous transition for the atmospheric turbulence realization temporally will require another method. The X _i ’s can be modified from being just random numbers to a continuous function of time for each mode. Thus, Equation (30) can be rewritten as,

$$\text{Wavefront}\left(\rho ,\theta \right)={\displaystyle \sum _{i=1}^M{X}_i{a}_i{K}_i\left(\rho ,\theta \right)}$$ E31

where the X _i (t) function here is generated by, first, creating a vector with a few random numbers with the zero mean unitary Gaussian distribution, as in Fig. 6.

Figure 6.

Sample vector with a few random numbers with the zero mean unitary Gaussian distribution.

Next, a spline curve is fit to this vector of a few random numbers, shown in Fig. 7, and this spline curve is now the X _i (t) temporal function for generating the wavefronts in the atmospheric turbulence simulation. Without this splining technique, the change between phase screens would be discontinuous and would not provide an accurate representation of the atmosphere for testing an adaptive optics system. In reality, the Earth’s atmosphere is a continuous medium. With this technique, the temporal transition of the wavefronts in the atmospheric turbulence simulation is continuous and smooth. Also, in conjunction with the use of Karhunen-Loeve polynomials, a statistically independent realization of the atmosphere is preserved.

It has been shown in various experiments that the first order aberrations, ie tip and tilt, are larger in magnitude and vary less with respect to time (Born & Wolf, 1997, Wilcox, 2005). To futher validate the Spline technique, one can take this into account by using a vector of fewer numbers than for the higher order aberrations for tip and tilt and the larger magnitude is taken care of by the Zernike-Kolmogorov residual errors used with the a _i ’s in Equation (31). Fig. 8 illustrates the temporal difference between the transitions of tip and tilt and those of some higher order aberrations. In the next section, a comparison between this technique and the Frozen Seeing model will be analyzed and discussed.

Figure 7.

Sample X _i (t) temporal function generated from a vector of a few random elements.

Figure 8.

Temporal and magnitude difference between (a) tip and tilt and (b) higher order aberrations.

4. Comparison of spline technique to the Frozen Seeing model

The method of generating atmospheric turbulence with temporal evolution as described in the previous section proposes various advantages compared to the Frozen Seeing model. The computational time to generate a phase screen of atmosphere of size NxN increases exponentially. Fig. 9 illustrates the number of seconds required to generate a phase screen of atmosphere using the Frozen Seeing model.

Figure 9.

Seconds to compute an NxN phase screen of atmosphere.

Figure 10.

Simulation of the Frozen Seeing model

Once the phase screen is generated, to be able to simulate the atmospheric turbulence on a SLM with the Frozen Seeing model, a subsection of that image of appropriate size is taken and used at the phase screen to represent the atmosphere at a moment in time. Then, that subsection is drifted across the large phase screen and that represents the next moment in time, to simulate the behavior of wind. This process is repeated until the edge of the NxN phase screen is reached, as shown in the illustration in Fig. 10.

One can clearly see that by generating atmospheric turbulence in this fasion will last for only a few seconds. Increasing the size of the large phase screen, N, would allow for a longer simulation, but the computational requirement to generate that phase screen would be a computational burden. In addition, with an NxN array, the number of bytes in that array will be N ². This will quickly lead to an image size of dozens of megapixels which will eventually lead to a software overflow. What can also be done is rather than drifting the subsection of the large phase screen across in a stright line is drifting in a circular motion about the large phase screen, but this will lead to a simulation of atmosphere that is very repetitious. Using the Spline technique outlined in the previous section, one can realistically simulate atmospheric turbulence for a longer period of time with far less computational requirements.

To compare the Frozen Seeing model to the Spline technique outlined in the previous section, each subsection of the larger phase screen can be analyzed with a single value decomposition (SVD) of the numerical values and calculate the Zernike coefficients, a _i ’s, of Equation (17) with M = 24. Fig. 11 illustrates the SVD (b) of a sample wavefront from a realization of atmosphere with a D/r ₀ = 2.25 (a), and the values of the SVD are listed in Table 3. Next, this process is repeated as the subsection is drifted across the large phase screen, to show the temporal transition of the a _i ’s.

Figure 11.

a) Sample wavefront of atmosphere with D/r ₀ = 2.25 and its (b) SVD of Zernike polynomials

The SVD representation in Fig. 11 (b) of the wavefront in Fig. 11 (a) has a fitted percent error of less than 2%. The a _i ’s progression over time are expected to change in a quasi-random fasion. It can be seen in Fig. 12 (a) and (b) that the temporal transitions of the a _i ’s resemble the temporal transitions as in the Spline technique, as shown in Fig. 13. Furthermore, the tilt components, a ₂ and a ₃ , are larger in magnitude than higher orders, which is consistant with the turbulence model outlined by the Zernike-Kolmogorov residual errors and the Spline technique.

i	ai	i	ai
1	-3.74173	13	0.344473
2	1.3423	14	-0.09171
3	0.384392	15	0.045092
4	-0.47419	16	0.179839
5	-1.51239	17	0.003466
6	0.398136	18	-0.01032
7	-0.2772	19	-0.14574
8	-0.4476	20	0.052195
9	-0.21973	21	0.121335
10	0.249914	22	0.354472
11	-0.06698	23	-0.18282
12	0.267201	24	-0.24076

Table 3.

Zernike polynomial coefficients that make up a sample representation of atmosphere with a D/r ₀ = 2.25

Figure 12.

The a _i ’s progression over time for (a) tip and tilt and (b) higher order Zernike terms in a simulation of atmospheric turbulence generated via the Frozen Seeing method.

Figure 13.

The a _i ’s progression over time for (a) tip and tilt and (b) higher order Zernike terms in a simulation of atmospheric turbulence generated via the Spline Technique.

By visual inspection, the two methods simulate atmospheric turbulence in a similar way. A statistical measure of the similarity of these two methods can be described by the cross-correlation of the respective a _i ’s.

Zernike order	Aberration	Average Cross-Correlation
a2	Tip	0.7004
a3	Tilt	0.7471
a4	Focus	0.6686
a5	Astigmatism X	0.7433
a6	Astigmatism Y	0.5937
a7	Coma X	0.5981
a8	Coma Y	0.6703
a9	Trefoil X	0.6909
a10	Trefoil Y	0.4878
a11	Spherical	0.5910
a12	Sec. Astigmatism X	0.5277

Table 4.

Average cross-correlation values for each a _i

After generating and analyzing ten realizations of atmosphere from the Frozen Seeing method, the average cross-correlation values are summarized in Table 4. An overall average of these cross-correlation values is 0.6381. This shows a consistancy between the generally accepted Frozen Seeing model and the new Spline technique outlined here.

5. System performance and results

The Holoeye LC2002 SLM device used in this example is a diffractive device that can directly modulate the phase of an incoming wavefront by π radians. In order to utilize the full 2π radian phase modulation on the impinging wavefront, one can set up a Fourier Filter and use either the +1 or -1 diffractive order through the rest of the system. The graphical user interface (GUI) of the software developed for controlling this system, written in Matlab, can be seen in Fig. 14.

Figure 14.

Graphical user interface for Holoeye Atmospheric Turbulence System

This software is capable of controlling any LC SLM. The Holoeye LC2002 is a device with 800x600 pixels and can accept a beam of 0.82“ in diameter. This software developed can set up the alignment of the diffraction orders and set alignment biases that may be entered to compensate misalignments in the optical components of the overall system for maximum performance. Different algorithms of generating turbulence can be used if desired and the parameters for the simulated telescope diameter and Fried parameter can control the severity of turbulence, as well. If desired, a secondary annular obscuration can be included in the simulation to simulate a telescope’s secondary mirror.

Figure 15.

Sample (a) and (c) wavefronts and (b) and (d) their corresponding PSFs due to atmospheric turbulence with a 0.4 meter telescope and an r ₀ of 1 cm.

Using the GUI developed, sample atmospheric conditions have been calculated and put on the SLM and then theire PSFs are measured with an imaging camera. The simulated atmospheric turbulence was calculated for a 0.4 meter telescope with seeing conditions having an r ₀ of 1 cm. Sample wavefronts from the simulation and their theoretical PSFs can be seen in Fig. 15.

The measured PSFs from the wavefronts in Fig. 15 (a) and (c) can be seen in Fig. 16 (a) and (b) and they are similar to that of the calculated PSFs in Fig. 15 (b) and (d), respectively. The 2-dimenstional cross-correlation factors between the two frames and their theoretical components are 0.9589 and 0.8638, respectively, showing that the system performs quite well and the measured and theoretical values are consistent with each other.

Figure 16.

PSF measurements of the two sample wavefronts in the optical system

Figure 17.

PSFs of (a) open-loop and (b) closed loop frame with x and y cross section plots.

At the Naval Research Laboratory, we have developed an AO system for use in astronomical applications (Restaino, S.R., et. al., 2008). We have simulated atmospheric turbulence with the system outlined in the previous sections and caused distortions on a laser beam for our AO system to correct. Simulating fairly reasonable seeing contitions with D/r ₀ = 1.5 lead to roughly a time-averaged Strehl ratio of 0.32. The Strehl ratio is a common way of measuring the effect that aberrations have on the imaging system (Born and Wolf, 1997). The typical definition of the Strehl ratio is the ratio of the peak intensity between the unaberrated system PSF and the system PSF with aberrations. Thus, a diffraction limited system, or a system limited only by the diffraction at the edge of the entrance pupil, will have a Strehl ratio of 1, and any aberration present in the system will cause the Strehl ratio to be less than 1.

Fig. 17 (a) and (b) show the PSFs of a frame taken during open-loop and closed-loop operation with their respective x and y cross sections. There is noticeable increase in peak intensity of the PSF and the other feature that is very important is the formation of the first ring of the Airy function in the corrected PSF. After closing the loop and allowing the AO system to begin correction, the time-averaged Strehl ratio for the simulation was increased to 0.84.

6. Summary

The method of generating atmospheric turbulence via the Spline technique is virtually the same as the Frozen Seeing method with the added feature of being far less computationally intensive on a computer system. This advantage can be exploited in the development of a software package that can drive any SLM to simulate atmospheric turbulence in almost any wavelength for any telescope diameter and adaptive optical and laser communication systems can be tested for performance evaluations. At the Naval Research Laboratory, a current system is being used with software written in the programming language Matlab and various tests are ongoing. Currently, two SLMs from Holoeye and Boulder Non-Linear Systems are being investigated and various wavelengths are being ustilized for different applications. Future work will include the investigation of other new liquid crystal devices as the field of liquid crystal technology is a very rapidly moving and growing field.