Atmospheric Propagation of High Energy Lasers: Thermal Atmospheric Propagation of High Energy Lasers: Thermal Blooming Simulation Blooming Simulation

. High Energy Laser (HEL) propagation through turbulent atmosphere is examined via numerical simulation. The beam propagation is modeled with the paraxial equation, which in turn is written as a system of equations for a quantum ﬂuid, via the Madelung transform. A ﬁnite volume solver is applied to the quantum ﬂuid equations, which supports sharp gradients in beam intensity. The atmosphere is modeled via a coupled advection-diffusion equation whose initial data have Kolmogorov spectrum. In this model the combined effects of thermal blooming, beam slewing, and deep turbulence are simulated.


INTRODUCTION
This article considers atmospheric effects on the propagation of continuous wave High Energy Lasers (HEL).The understanding of laser propagation is important in an ever increasing array of applications including targeting, wireless communication, remote sensing, the creation of Bose-Einstein Condensate, gravity-wave measurements, and many more [1,2,3,4].For applications involving atmospheric propagation, continuous wave lasers exist with powers in the kilowatts [5], a number which is certain to increase with time.In laboratory settings, pulse lasers have been developed with powers reaching the staggering petawatts [6,7].In this work we develop a numerical solver for the propagation of continuous wave HEL, which is capable of handling three important environmental interactions: scintillation by random media, thermal blooming, and beam slewing (by wind).In the context of this solver the combined effect of these interactions is presented.
It is unreasonable to directly simulate the dynamics of the atmosphere at the scale of laser wavelengths.Instead, as is typical, we use a statistical description of the atmosphere.To model a turbulent atmosphere, we use a random temperature field with prescribed spectrum.This temperature field then induces a random refractive index.In this work, the beam propagates through a random refractive index for its entire length.This contrasts much of the previous work on propagation through random media which focuses discretely placed phase screens [8,9,10,11,12,13].The statistics of the temperature field are designed to match an expression for the expected intensity, as formulated by [14].
Even in a statistical description of the atmosphere, one may take into account the effect of the beam heating the air it passes through.The resulting change in refractive index feeds back upon the beam, an effect known as thermal blooming.Thermal blooming is a well documented phenomenon [15], whose relevance to laser propagation in the atmosphere gains importance as laser's become more powerful.As a simple coupling, we append an advection-diffusion equation for the temperature, in which heating from the laser appears as a forcing term.The temperature feeds back on the beam propagation via the refractive index.The beam's propagation is modeled using the paraxial, or Schrödinger, equation.This standard model can be derived directly from Maxwell's equations when backscatter effects are neglected [16].The numerical method used here is based in reformulation of the paraxial equation, via the Madelung transformation [17,18,19,20].This formulation divides the complex field into two real valued fields that are smoother than the original complex field.The transformed equations are similar to the equations of isentropic gas dynamics, but have an extra term called the quantum pressure.This approach was previously used in chemistry [21] and has been suggested for optics [18].The Madelung transformed fluid equations are evolution equations in propagation distance, whereas the temperature equation is an evolution equation in time.We evolve these equations alternatively in space-time, using a simple splitting scheme.
The remainder of the paper is organized as follows.In section 2, the equations modeling the beam propagation and temperature field are presented.In section 3, the numerical methods used for simulation are presented.Also in this section is a description of the technique used to generate the random smooth fields and an example of a beam's propagation through random media with thermal blooming and a side-wind.In section 4, we conclude and present future research areas.

FORMULATION
In this section, we present the model equations used to study High Energy Laser (HEL) propagation through the air.These equations are the paraxial equation for the beam propagation, coupled with a forced advection-diffusion equation for the temperature fluctuations in the atmosphere.The forcing term in the temperature equation models the heating of the beam, and allows for the study of thermal blooming.The advection term models allows for the modeling of side-wind.Random initial temperature distributions serve as the model for a turbulent atmosphere.
We choose to approximate the beam's propagation by the paraxial equation, as in [20].As an approximation to Maxwell's equations this amounts to neglecting backscattering and polarization, and assuming small changes in the refractive index.From this point forward, consider x ∈ R 2 to be the x, y-coordinate directions parallel to the plane containing the aperture and z ∈ [0, Z] to be the coordinate in the propagation direction.We write the refractive index as, ( 1) where n 0 is the average value of the refractive index and n 1 a small real valued perturbation, (2) We expect the change in the refractive index due to thermal blooming to be about 10 −9 [9].
The paraxial equation for the beam propagation is where H is the Laplacian in the xy-plane and P(x, z) is the potential, defined as ( 4) This equation has been numerically simulated previously using, for example, the Fourier-split step method [8,9,10,13].For propagation through turbulent media, the field V often develops sharp gradients, for which pseudospectral methods are ill-suited.We handle these gradients by working in the Madelung transformed paraxial equation [21].
From the field, V , in equation ( 3), the beam's intensity, I, is defined as, ( 5) We will denote ρ(x, z) as the rescaled intensity.
Equation ( 3) can be rewritten in a way to mimic a fluid equation.We will rewrite the field where ρ is the rescaled intensity and θ is the phase.Both ρ and θ are real valued fields.
Madelung introduced this transform in order to move from a linear (and nonlinear) Schrödinger equation to hydrodynamic type system of equations [22,23].Note that absorption is still included in the simulation due to equation (5).The Madelung transform of equation ( 3) yields a fluid-like form of the paraxial equation, in which v = ∇ H θ .In equation (7), Im(n 0 n 1 ) has been neglected.Since θ is multivalued, the gradient is numerically evaluated as, An alternative way would be to follow the method in [24].Note that equation (7a) is invariant to rescaling of ρ, therefore for potentials which are independent of ρ, a low power laser will propagate along the same path as a high power laser.When thermal blooming is included, the potential will depend on ρ and this is no longer the case.
As the power of the laser increases, with all other aspects fixed, the amount of energy the laser loses to the medium in the form of heat also increases.This heating feeds back into the dynamics of the laser, at leading order, via the effect of temperature on the refractive index.Since the speed of light is much larger than the speeds associated with the change in temperature, the equation system ( 7) can be solved using a fixed temperature.The solution of (7) will then play the role of a known forcing in an evolution equation for temperature.We choose to model the temperature fluctuation T with a convection-diffusion equation, ( 9) where u(x, z,t) is some known velocity profile, modeling wind, where ε is the absorptivity, I the power of the beam per unit area [9].Such a model is reasonable in the case where the timescale over which the laser operates is small relative to the time it takes for convective-driven fluid motion and that such motions are small relative to the background wind (otherwise one should also simulate the velocity field).
Note that the absorptivity is related to the complex part of the refractive index by the Beer-Lambert Law, [25], (10) The above equation is used to determine the value of Im(n 0 ).To complete the connection between equation ( 7) and ( 9) one needs to specify a closure between refractive index and temperature.We posit that the variation of the refractive index depends linearly on the temperature fluctuation, as in [9], (11) n 1 = −10 −6 T.

NUMERICAL SIMULATION OF HEL THROUGH RANDOM MEDIA
In order to numerically simulate thermal blooming, we need to solve two partial differential equations.The two partial differential equations are of different type and will be simulated with different numerical algorithms.The beam propagation will be simulated by a hyperbolic finite volume solver, see section 3.1, whereas the advection-diffusion equation will be simulated using an upwind finite difference solver.Since the beam propagation happens on a much faster time scale than the heating of the medium, it is possible to separate the two system of equations.
The refractive index is frozen and the light is assumed to travel instantaneously from the source to the target.
3.1.Numerical Method.Equations (7a) and (7b) will be solved using the package CLAW-PACK [26].CLAWPACK requires an approximate Riemann solver for equations (7).Since the left hand side of equation (7b) is a rescaled inviscid Burgers equation, The solver is straight forward to implement.
For the transonic rarefaction wave solution, v r < 0 and v l > 0, we will use the entropy fix described in [27], where A ± ∆v are the fluxes of the velocity in the positive and negative direction and A ± ∆ρ are the fluxes of the density in the positive and negative direction.The extension to two dimensions is straightforward.The source term is approximated by the five-point stencil finite difference method.To prevent the numerical differentiation of the quantum pressure from being unbounded the following substitution will be used, where √ ρ = e R .This prevents a numerical small divisor problem.
To numerically simulate the temperature equation ( 9), we use a finite difference method.The advection term is calculated using first order upwind and the diffusion by the five-point stencil finite difference method.The time stepping algorithm used is variable order Adams-Bashforth-Moulton PECE, see [28], by the command ode113 in MATLAB.Note that the electromagnetic field and the the temperature field can be solved using different grids.The fields are interpolated using linear interpolation.

Generating refractive index for turbulent air.
The traditional way to numerically simulate propagation through turbulent air is to use discretely placed phase screens, see [29].While phase screen method produces good results, it has not been rigorously shown that it gives the same intensity profile as propagation through a continuous random medium.We instead generate a random refractive index with Kolmogorov statistics.The change in refractive index for turbulent air is modeled as a passive scalar.Passive scalars have a similar correlation as Kolmogorov-type turbulence, see [30].We generate the refractive index by following [31], a method for generating fractional Brownian motion, (1) Generate a field of random complex numbers in Fourier space, by setting each Fourier coefficient to, where U k 1 ,k 2 ,k 3 and V k 1 ,k 2 ,k 3 are independent standard normal random variables for all k 1 , k 2 , k 3 .Since we want to find a real valued field, we use the complex conjugate to find Fourier coefficients in the variables x, y and z respectively.(2) Scale each Fourier coefficient, so that the three dimensional spectrum is set to, see [32], 3 the above parameters are • L, integral length scale, here 50m, • l K , Kolmogorov micro-scale, here 1mm, • c L , a positive constant, here 6.78, • p 0 , a positive constant, determining the power of the three dimensional energy function spectrum as κ → 0, p 0 = 4, • β = 5.2 and c η = 0 are positive constants.
(3) Rescale f from equation (14), in order to obtain the desired energy function spectrum.(16) f (4) Use the inverse fast Fourier transform to obtain the perturbation of the refractive index in real space: This method is described in more detail for fractional Brownian noise in [33].This construction creates a refractive index correction n 1 (x, y, z) which is periodic, discrete in space and Gaussian.
In figure 1, a slice of one such realization is depicted in the xy-plane.The average has been deducted from the field to highlight the perturbations.Even though the horizontal gradient is small, the effect on propagation is large.
To test the randomness of a particular realization, we use the correlation function, F(r) = n 1 (x, y, z)n 1 (x, y, z + r)dxdydz n 1 (x, y, z)n 1 (x, y, z)dxdydz , which is approximated by, ( 19) The approximate correlation function is shown in the left hand side of figure 1.For large distances, the approximate correlation function will oscillate around zero.If more resolution is used in z-direction, the size of the oscillations will decrease.

Numerical Results.
As discussed above, our numerical method alternately solves (7) with a finite volume method and (9) with a finite difference method.The numerical parameters are M x = M y = 128 equally spaced points in the transverse directions, M z = 61 z-planes, with ∆z=82m and ∆t=0.007s.The intensity is reported in units of 100 kW/m 2 .In the top row, from left to right, t = 0s, 0.015s, and 0.036s; in the bottom row t = 0.043s, 0.05s and 0.65s.
The simulated physical parameter regime is similar to that of [9] In figure 2, we show the results of laser propagating through 5 km of air with thermal blooming.The initial temperature field was generated using the method described in section 3.2.
The temperature field evolves according to equation (9).At time t = 0, the laser propagates through a random field and generates hotspots of very high intensity, e.g.scintillation.As time progresses, the intensity field smooths due thermal blooming -the hotter the spots the more pronounced the thermal blooming.For large time, figure 2, thermal blooming is the dominating effect on propagation.In the last time slice, the intensity distribution retains some asymmetry due to a combination of the initial temperature field and the sidewind.

CONCLUSIONS
HEL propagation was simulated using the paraxial equation coupled with an advectiondiffusion equation for the temperature field.The quantum fluid formulation, derived by [23], was employed.This fluid formulation enables the use of shock capturing finite volume methods, allowing for simulation of steeper gradients than traditional solver methods, for example split step Fourier method.The effects of beam slewing by wind, scintillation by random media, and thermal blooming were simulated in concert.We observe that thermal blooming helps to mitigate the scintillation of the random medium, and beam slewing counteracts thermal blooming.
Disclaimer: This report was prepared as an account of work sponsored by an agency of the United States Government.Neither the United States Government nor any agency thereof, nor any of their employees, make any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights.Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof.The views and opinions of authors expressed herein do not necessarily state or reflect those of the

FIGURE 1 .
FIGURE 1. LEFT: The correlation function in z direction, of the refractive index, is depicted as a function of distance r, in units of the integral lengthscale L. RIGHT: A realization of the random perturbation to the refractive index, n 1 , for Kolmogorov turbulence is depicted for at one point in z.The average value of the refractive index in the plane has been subtracted.The scale is in 10 −4 .

FIGURE 2 .
FIGURE 2. Intensity of the beam after propagation through random temperature field with sidewind v x = −0.5m/sand including the effect of thermal blooming.