Robert M. Manning

Glenn Research Center, Cleveland, Ohio

High Energy Laser Beam Propagation in the 

Atmosphere: The Integral Invariants of the 

Nonlinear Parabolic Equation and the 

Method of Moments

NASA/TM—2012-217634

July 2012
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientific and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Officer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of significant 

scientific and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443–757–5803

• Phone the NASA STI Information Desk at

443–757–5802

• Write to:

 STI Information Desk

 NASA Center for AeroSpace Information

 7115 Standard Drive

 Hanover, MD 21076–1320
Robert M. Manning

Glenn Research Center, Cleveland, Ohio

High Energy Laser Beam Propagation in the 

Atmosphere: The Integral Invariants of the 

Nonlinear Parabolic Equation and the 

Method of Moments

NASA/TM—2012-217634

July 2012

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135
Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://www.sti.nasa.gov

Level of Review: This material has been technically reviewed by technical management. 

This report contains preliminary findings, 

subject to revision as analysis proceeds.
NASA/TM—2012-217634 1

High Energy Laser Beam Propagation in the Atmosphere: 

The Integral Invariants of the Nonlinear Parabolic 

Equation and the Method of Moments 

Robert M. Manning 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Abstract 

The method of moments is used to define and derive expressions for laser beam deflection and beam 

radius broadening for high-energy propagation through the Earth’s atmosphere. These expressions are 

augmented with the integral invariants of the corresponding nonlinear parabolic equation that describes 

the electric field of high-energy laser beam to propagation to yield universal equations for the 

aforementioned quantities; the beam deflection is a linear function of the propagation distance whereas 

the beam broadening is a quadratic function of distance. The coefficients of these expressions are then 

derived from a thin screen approximation solution of the nonlinear parabolic equation to give 

corresponding analytical expressions for a target located outside the Earth’s atmospheric layer. These 

equations, which are graphically presented for a host of propagation scenarios, as well as the thin screen 

model, are easily amenable to the phase expansions of the wave front for the specification and design of 

adaptive optics algorithms to correct for the inherent phase aberrations. This work finds application in, for 

example, the analysis of beamed energy propulsion for space-based vehicles. 

1.0 Introduction 

When delivering large amounts of power through the Earth’s atmosphere via millimeter or infrared 

beams (i.e., laser beams or beams formed at the output of a millimeter wave antenna system), many 

propagation mechanisms must be addressed that may be potentially deleterious to such power 

transmission. The most obvious one is the ever-present random variation of the atmospheric refractive 

index due to local temperature variations known as ‘turbulence’. This naturally occurring phenomena is 

driven by thermal convection of heat from the Earth’s surface; once the resulting air motion exceeds a 

critical value of velocity, laminar flow essentially evolves into turbulent flow and fluctuations in the 

temperature distribution becomes statistically random (Ref. 1). These temperature fluctuations then act 

directly on the prevailing refractive index, thus rendering the refractive index a random quantity. These 

refractive index variations randomly focus and defocus the intervening electromagnetic wave field. Thus, 

the atmosphere can be considered to be composed of ‘lenses’ of random focusing and defocusing 

characteristics that, due to the gross atmospheric motion due to wind, move across the beam. This gives 

rise to many beam quality variations; the major ones being beam broadening and beam steering. The 

statistical analysis and modeling of this type of atmospheric propagation as a long and rich history and 

has resulted in analytical descriptions for the impact of turbulence on the operation of systems relying on 

such beam propagation. Many models and descriptions exist for the ‘engineering analysis’ of the 

operation of transmission systems that rely on the propagation of electromagnetic beam propagation in the 

atmosphere (Ref. 2). 

The scenario discussed above may be considered as ‘passive’ electromagnetic wave propagation, i.e., 

the wave field moves through an atmosphere the refractive index of which is determined by other sources, 

not the field itself. However, as the energy density of the beam increases, absorption of the beam energy 

by atmospheric gas components results in local heating of the atmospheric which does indeed act directly 

on the refractive index causing it to decrease in value. This thermal change of the refractive index field 
NASA/TM—2012-217634 2

then acts on the electromagnetic wave field causing it to also change, and so on. The propagation scenario 

now becomes an ‘active’ one, whereby the propagating field modifies the very medium it which it exists. 

This heating process is called ‘thermal blooming’ and substantially differs from that of the passive 

propagation discussed earlier. Here, a ‘thermal lens’ is created within the atmosphere by the heating due 

to the energy density of the beam. This ‘self-action’ of the beam will not only bend the beam into regions 

of higher refractive index (beam steering), but convection within the atmospheric fluid will also arise 

which is the source of self induced turbulent flow of the medium. The situation is further complicated 

when one includes the effects of atmospheric wind and aerosols and the abovementioned passive 

propagation effects. Defocusing and other such associated nonlinear thermal blooming distortions of the 

beam cross-section will then result. In extreme cases of very large energy densities, the Kerr effect will 

arise and the propagating beam can subsequently break up into smaller beams, or filaments, which 

severely constrains the amount of energy density that the beam will be able to possess as it travels through 

the atmosphere. Unlike the situation of passive propagation, the thermal blooming mechanism introduces 

nonlinearities into the analysis of the phenomena that substantially complicates a complete mathematical 

description. Complete analyses of these types of propagation scenarios can only be done numerically, 

which was a major activity within the United States and Russia in the late 1980s. Other than the usual 

‘order-of-magnitude’ estimates using the equations of fluid mechanics and wave propagation, only 

numerical modeling of the effects of atmospheric thermal nonlinearities abound in the literature. 

Analytical treatments appropriate for an engineering analysis and assessment of atmospheric propagation 

systems encountering thermal blooming have been lacking and those that do exist hold only for the 

geometrical optics case (Refs. 3 to 5). 

Two of the more important performance parameters of a high-energy laser-beaming scenario are the 

accuracy to aim the beam in the desired direction and the size of the beam at the intended target. In most 

cases, the beam will be focused at the target so the relevant quantity is the size of the focal spot on the 

target. It is the purpose of this memorandum to derive expressions for the beam deflection from the 

intended axis and the broadening of the beam focal spot at the target using only the information that can 

be gleaned for the structure of the nonlinear parabolic equation that describes the propagation process in 

the presence of thermal nonlinearities without actually solving the equation in its entirety. This will be 

accomplished by invoking the integral invariants that exist for the nonlinear parabolic equation as well as 

employing the method of moments to express the statistical moments of the electric field that are used to 

define the deflection as well as the radius of the laser beam. Once this has been done, a thin screen 

approximation solution will be given for the propagating wave field of the laser beam through the Earths 

atmosphere and these results will be applied to derive analytical and graphical results for the beam 

parameters focused to a target above the atmospheric layer. This thin screen model will allow for the 

design and specification for the adaptive correction of the phase aberrations that prevail in such 

propagation situations, and the expressions for the deflection and radius will allow for the evaluation of 

their corrected values. 

2.0 The Moments of the Laser Beam Energy Density and the Nonlinear 

Parabolic Equation 

Given the energy density E x( ) , y,z

2

 of the electromagnetic field distribution of a laser beam in a 

plane 

r⊥ = xxˆ + yyˆ transverse to the direction of propagation at a distance z from the origin, one can 

define the mnth moment of the deposition of this quantity in the transverse plane, viz, 

γ mn ( )z = xm yn E x( ) , y,z

2

−∞

∞

∫ −∞

∞

∫ dxdy (1) 
NASA/TM—2012-217634 3

Thus, γ00(z) is the total beam energy at the distance z, γ10(z) is the first moment or centroid of the beam 

energy from the z-axis in the x-direction, γ01(z) as that in the y-direction, γ20(z) as the second moment of 

the beam energy in the x-direction, etc. Hence, using these spatial moments of the energy density in a 

plane, one can now form the vector quantity 



Δ⊥ ( )z ≡ 1

γ 00 ( )z

γ10 ( )z xˆ + γ 01 ( ) ( )z yˆ (2) 

that can be defined as the normalized overall transverse displacement of the beam energy from the z-axis 

(intended direction of propagation) at the distance z, and 

aeff

2 ( )z ≡ 1

γ 00 ( )z

γ 20 ( )z + γ 02 ( ) ( )z (3) 

that can be defined as the square radius of the spread of the beam energy about the z-axis. This particular 

moment can be centralized by forming the quantity 

reff

2 ( )z = aeff2 ( )z − 

Δ⊥ ( )z (4) 

which now gives the radius of the beam energy (i.e., the beam radius) about the axis shifted by 

Δ⊥ ( )z

from the intended direction. 

The analysis of these quantities can be a rather difficult problem in the full nonlinear case of high 

energy laser beam propagation in an absorbing medium which gives rise to a host of thermal 

nonlinearities collectively known as thermal blooming. It thus becomes of interest to ascertain just how 

much of the behavior of these performance parameters can be gleaned from the fact that the equation 

governing the electric field of the laser beam in an absorptive nonlinear medium is given by 

2ik ∂E

∂z

+ ∇⊥

2 E + k2

εT T E 2

( ) E = 0 (5) 

without fully solving specific problem scenarios. Here, ∇⊥

2 is the Laplacian taken along the direction 

r

transverse to the direction of propagation z. Equation (5) is the well-known parabolic equation applied to 

a propagation medium whose variation of permittivity Δε is governed by the change in temperature T of 

the medium from its nominal value which is determined by the energy density E 2 of the field, i.e., 

Δε T, E 2

( ) = εT T E 2( ) (6) 

where εT ≡ ∂ε ∂T is the variation of the permittivity with respect to T; in the case of the Earth 

atmosphere, εT < 0 . The functional dependency of the temperature variation on the field energy density 

is given by the application of the conservation of energy to the atmospheric heat budget, given in its 

entirety by 

ρCp

∂T

∂t

+ V ⋅



∇T ⎛

⎝

⎜ ⎞

⎠

⎟ = κ∇2

T + α

c

8π

E 2 (7) 

where ρ is the density of the atmosphere, Cp is its specific heat at constant pressure, κ is its thermal 

conductivity, 

V is the atmospheric wind velocity, and α is the absorption coefficient. (Strictly speaking, 
NASA/TM—2012-217634 4

one should also include the set of Navier-Stokes equations in the Boussinesq approximation to account 

for the viscous and Archimedean forces that occur during the thermally induced motion of the 

atmosphere. However, for the purposes of this discussion and subsequent analysis, a description at this 

level is not required.) The implications of the parabolic equation of Equation (5) applied to this nonlinear 

problem will now be investigated for the evolution of the moments of the beam parameters discussed 

above. 

3.0 The Invariants of the Nonlinear Parabolic Equation 

Two integral invariants (sometimes called the ‘integrals of motion’) of the differential equation of 

Equation (5) can easily be derived (Refs. 6 to 8). First, multiplying Equation (5) by E*

 and subtracting 

from it the complex conjugate of the same equation gives 

2ik E∗ ∂E

∂z

+ E ∂E∗

∂z

⎛

⎝

⎜

⎞

⎠

⎟ + E∗

∇⊥

2 E − E∇⊥2 E∗ = 0 (8) 

Using the fact that E 2 ≡ EE∗ and simplifying where possible yields 

2ik ∂ E 2

∂z

+ ∇⊥ ⋅ E∗ 

∇⊥E − E∇⊥E∗ ( ) = 0 (9) 

Finally, integrating Equation (9) in the transverse plane yields and applying Gauss’s Law to the second 

member of the relationship and letting the surface approach infinity, at which the surface terms approach 

zero, one finally obtains 

E 2

−∞

∞

∫ −∞

∞

∫ dxdy = const = γ 00 (10) 

which is just a statement of the conservation of energy for any transverse plane. (Absorption of the wave 

field is present which, after all, accounts for the heating and subsequent nonlinear action of the medium. 

However, the action of the absorption on the overall energy density of the laser beam itself has been 

neglected in Equation (5). Incorporating it into the subsequent development is straightforward but will not 

be done here.) 

Second, multiplying Equation (5) by ∂E∗ ∂z and adding to it the complex conjugate of the same 

equation gives 

∂E∗

∂z

∇⊥

2 E +

∂E

∂z

∇⊥

2 E∗ + k2

εT T E 2

( ) ∂ E 2

∂z = 0 (11) 

As done above, one now integrates this relationship over the transverse plane. Each of the first two terms 

on the left side of Equation (11) can then be simplified by first using Green’s first identity and letting the 

surface terms vanish for the same reason as employed above: 

− ∇⊥

∂E∗

∂z ⋅



∇⊥E + ∇⊥

∂E

∂z

⋅



∇⊥E∗ ⎛

⎝

⎜

⎞

⎠

⎟ dx dy

−∞

∞

∫ −∞

∞

∫ + k2

εT T E 2

( ) ∂ E 2

∂z −∞

∞

∫ −∞

∞

∫ dxdy = 0 (12) 

Defining the quantity G E 2

( ) such that 
NASA/TM—2012-217634 5

T E 2

( ) =

∂G E 2

( )

∂ E 2 (13) 

and substituting it into Equation (12) and simplifying yields 

− ∂

∂z



∇⊥E ⋅



∇⊥E∗ ( ) dx dy

−∞

∞

∫ −∞

∞

∫ + k2

εT

∂

∂z

G E ∗

( ) −∞

∞

∫ −∞

∞

∫ dxdy = 0 (14) 

which becomes 



∇⊥E ⋅



∇⊥E∗ − k2εTG E 2

( ) ⎡

⎣ ⎤⎦

dx dy

−∞

∞

∫ −∞

∞

∫ = const ≡ H (15) 

which is essentially the Hamiltonian for the propagating wave field. 

4.0 Conservation Laws for the Nonlinear Parabolic Equation 

The two integral invariants given above should now be augmented with the conservation laws that 

prevail quantities that issue from Equation (5). A streamlined method to accomplish this is to employ 

field theory and find the Lagrangian corresponding to Equation (5) and from it obtain the energy 

momentum tensor from which the conservation quantities and laws would follow (Ref. 9). However, a 

direct approach will be used using Equation (5) directly. To this end, the first such conservation law was 

already derived and is given by Equation (9). Defining 



S⊥ ≡ − i

2k

E∗ 

∇⊥E − E



∇⊥E∗ ( ) (16) 

which is the energy flux vector of the field, one has 

∂ E 2

∂z

+ ∇

ρ ⋅



S⊥ = 0 (17) 

This is one of two prevailing conservation laws. 

Consider now the other equation which governs the energy flux vector 

S⊥ ; the idea here is to obtain 

an equation for 

S⊥ analogous to that of Equation (17) for E 2 . Using Equation (16), 

∂



S⊥

∂z = − i

2k

⎛

⎝

⎜ ⎞

⎠

⎟ ∂E∗

∂z



∇⊥E + E∗ ∇⊥

∂E

∂z

⎛

⎝

⎜ ⎞

⎠

⎟ − ∂E

∂z



∇⊥E∗ − E



∇⊥

∂E∗

∂z

⎛

⎝

⎜

⎞

⎠

⎟

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

 (18) 

Employing Equation (5) and simplifying gives 

∂



S⊥

∂z = 1

2k

⎛

⎝

⎜ ⎞

⎠

⎟

2

−K∗ 

F + E∗ 

∇⊥K − K 

F∗ + E



∇⊥K∗ ⎡

⎣ ⎤⎦ (19) 

where 
NASA/TM—2012-217634 6



F ≡ ∇⊥E, K ≡ 

∇⊥ ⋅



F + k2εT TE (20) 

Making use of the identity 



∇⊥ E∗ ( ) K = 

F∗K + E∗ 

∇⊥K (21) 

as well as its complex conjugate allows Equation (19) to become 

∂



S⊥

∂z = 1

2k

⎛

⎝

⎜ ⎞

⎠

⎟

2

−2



FK∗ − 2



F∗K + ∇⊥ E∗ ( ) K + ∇⊥ EK∗ ( ) ⎡

⎣ ⎤⎦ (22) 

Two more identities are now needed: 

∇⊥

2 E 2 = 2



F ⋅



F + E∗ 

∇⊥ ⋅



F + E



∇⊥ ⋅



F∗ (23) 

and 



F



∇⊥F∗ + 

F∗ 

( ) ∇⊥F

i

= ∂

∂xj

FiFj

∗ ( ) +

∂

∂xj

Fi

∗ ( ) Fj − ∂

∂xi



F ⋅



F∗ ( ) (24) 

(The latter comes from the fact that 

F × ∇⊥ × 

F∗ = 0 since 

F = ∇⊥E , etc.) After a bit of algebraic 

manipulation, one finds that 

∂



S⊥

∂z

⎛

⎝

⎜

⎞

⎠

⎟

i

= 1

2k

⎛

⎝

⎜ ⎞

⎠

⎟

2

−2 ∂

∂xj

FiFj

∗ + Fi∗ ( ) Fj +

∂

∂xj

∇⊥

2 E 2

( )δij

⎡

⎣

⎢

⎢

+2k2

εT

∂

∂xj

T E 2

( ) − T ∂

∂xj

E 2 ⎛

⎝

⎜

⎞

⎠

⎟ δij

⎤

⎦

⎥

⎥

 (25) 

Finally, using Equation (13) and simplifying, Equation (25) gives 

∂



S⊥

∂z

⎛

⎝

⎜

⎞

⎠

⎟

i

= 1

k

⎛

⎝

⎜ ⎞

⎠

⎟

2 ∂

∂xj

1

4

∇⊥

2 E 2

( )δij − 1

2

⎡

⎣

⎢



( ) ∇⊥E

i



∇⊥E∗ ( )j

− 1

2



∇⊥E∗ ( )i



( ) ∇⊥E

j

+

1

2

k2εT T E 2 − G E 2

( ) ( ) δij

⎤

⎦

⎥

 (26) 

Equation (26) can be written as the tensor relation 

∂



S⊥

∂z − 

∇⊥ ⋅



Q = 0 (27) 

where 

Q is the stress tensor 
NASA/TM—2012-217634 7



( ) Q ij ≡ 1

k

⎛

⎝

⎜ ⎞

⎠

⎟

2 1

4

∇⊥

2 E 2

( )δij − 1

2

⎡

⎣

⎢



( ) ∇⊥E

i



∇⊥E∗ ( )j

− 1

2



∇⊥E∗ ( )i



( ) ∇⊥E

j

+

1

2

k2εT T E 2 − G E 2

( ) ( ) δij

⎤

⎦

⎥

 (28) 

The conservation laws of Equation (17) and (27), as well as Equations (10) and (15), which prevail 

for Equation (5) are now sufficient to obtain expressions for the beam parameters specified by Equations 

(2) to (4). 

5.0 The Displacement and Broadening of a Laser Beam Experiencing 

Thermal Blooming 

The displacement from bore sight of a laser beam which induces thermal nonlinearities during 

propagation can now be easily calculated using Equation (2) involving the first order spatial moment of 

the energy density. To this end, differentiating Equation (2) (noting the constancy of γ00 by Equation 

(10)), and using Equation (17) yields 

d



Δ⊥ ( )z

dz = − 1

γ 00

x ∇⊥ ⋅

 ( ) S⊥ xˆ + y 

∇⊥ ⋅

 ( ) S⊥ yˆ ( ) −∞

∞

∫ dx dy

−∞

∞

∫ (29) 

Integrating by parts and noting that the surface terms approach zero at infinity gives 

d



Δ⊥ ( )z

dz = 1

γ 00



S⊥

−∞

∞

∫ dx dy

−∞

∞

∫ ≡ 

A z( ) (30) 

This is sufficient to describe the displacement in that 

d2 

Δ⊥ ( )z

dz2 = 1

γ 00

∂



S⊥

∂z −∞

∞

∫ −∞

∞

∫ dxdy = 1

γ 00



∇⊥

−∞

∞

∫ −∞

∞

∫ ⋅



Qdxdy = 0 (31) 

upon using Equation (27); the last term vanishes upon using Gauss’s Theorem and, once again, letting the 

fields and their derivatives vanish at infinity. Thus, 

A z( ) = 

A = const and the overall beam displacement 

is given simply by integrating Equation (30), i.e., 



Δ⊥ ( )z = Az (32) 

where it is assumed from hear on that 

Δ⊥ ( ) 0 = 0 , i.e., the beam is initially pointed down its intended 

axis. Thus, the beam deflection from its intended axis is a linear function of distance travelled. 

Similarly, the associated overall beam broadening reff

2 can be found by first differentiating 

Equation (3) and using Equation (17), 

daeff

2 ( )z

dz = − 1

γ 00

x

2 + y2 ( ) −∞

∞

∫ 

∇⊥ ⋅

 ( ) S⊥ dx dy

−∞

∞

∫ (33) 

Integrating by parts gives 
NASA/TM—2012-217634 8

daeff

2 ( )z

dz = 2

γ 00



r⊥ ⋅



S⊥

−∞

∞

∫ dx dy

−∞

∞

∫ ≡ B z( ) (34) 

Differentiating once again and using Equation (27) yields 

d2aeff

2 ( )z

dz2 = 2

γ 00



r⊥ ⋅ 

∇⊥ ⋅

 ( ) Q

−∞

∞

∫ dx dy

−∞

∞

∫ (35) 

Integrating by parts and letting the surface term vanish, 

d2aeff

2 ( )z

dz2 = − 2

γ 00

nˆ ⋅

 ( ) Q

−∞

∞

∫ dx dy

−∞

∞

∫ = − 2

γ 00

Qii

i=1

2

∑

−∞

∞

∫ dx dy

−∞

∞

∫ (36) 

where nˆ is the surface normal. In order to assure that this result is a constant and that no further 

derivatives are needed, one must evaluate the required trace of the tensor 

Q . Using Equation (28), one 

has 

Qii

i=1

2

∑ = 1

2

1

k

⎛

⎝

⎜ ⎞

⎠

⎟

2 

∇⊥ ⋅ E∗ ∇⊥E + E



∇⊥E∗ ( ) − 2



∇⊥E ⋅



∇⊥E∗ + 2k2εT T E 2 − G E 2

( ) ( ) ⎡

⎣⎢ ⎤

⎦⎥ (37) 

Introducing this into Equation (36) and using the Gauss Theorem to rid of the divergence term finally 

gives 

d2aeff

2 ( )z

dz2 = − 1

γ 00

⎛

⎝

⎜

⎞

⎠

⎟

1

k

⎛

⎝

⎜ ⎞

⎠

⎟

2

−2



∇⊥E ⋅



∇⊥E∗ + 2k2εT T E 2 − 2k2εTG E 2

( ) ⎡

⎣ ⎤⎦ −∞

∞

∫ dx dy

−∞

∞

∫ (38) 

This expression can be stated in terms of the invariant of Equation (15), viz, 

d2aeff

2 ( )z

dz2 = 1

γ 00

⎛

⎝

⎜

⎞

⎠

⎟

1

k

⎛

⎝

⎜ ⎞

⎠

⎟

2

2H − 2k2εT T E 2 − 2G E 2

( ) ⎡

⎣ ⎤⎦ −∞

∞

∫ dx dy

−∞

∞

∫ ⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

 (39) 

Hence, this expression is constant so long as the integrals vanish within the brackets. Thus, in the case of 

the Kerr effect, from Equation (6), one has that the temperature is replaced by T ⇒ A E 2 where A is 

some parameter independent of the electric field intensity, then by Equation (13), G E 2

( ) = ( ) 1 2 A E 4

and the integrand is identically equal to zero. However, in the case of thermally induced nonlinearities, 

the temperature field is related to the electric field intensity through Equation (7). In the case of forced 

convection due to wind and where large Peclet numbers prevail (see Appendix A), ∂T ∂x = β E x( ) , y,z

2

(for a wind velocity in the x-direction), and this will give rise to much more complicated functional 

dependencies. A general analysis of this problem is beyond the scope of this particular work but as shown 

in Appendix A, in the case of forced convection, 

T E 2 − 2G E 2

( ) ⎡

⎣ ⎤⎦ −∞

∞

∫ dx dy

−∞

∞

∫ ≈ 0 (40) 
NASA/TM—2012-217634 9

Hence, 

d2aeff

2 ( )z

dz2 = 1

γ 00

⎛

⎝

⎜

⎞

⎠

⎟

1

k

⎛

⎝

⎜ ⎞

⎠

⎟

2

2H = const ≡ C (41) 

Integrating this expression and remembering Equation (34) gives for the overall effective beam radius 

aeff

2 ( )z = 1

2

Cz2 + B( ) 0 z + aeff

2 ( ) 0 (42) 

and from Equations (4) and (32), for the effective radius about the deflected axis, 

reff

2 ( )z = 1

2

Cz2 + B( ) 0 z + aeff

2 ( ) 0 − 

( ) Az 2

 (43) 

Hence, Equations (32) and (43) give the universal behavior of the deflection and broadening of a laser 

beam governed by Equation (5) without solving the entire propagation problem. The constants 

A , B(0), 

and C are related to the beam parameters via Equations (30), (34), and (41), respectively. 

As a quick check on the validity of Equation (43), consider a Gaussian beam wave propagating in free 

space. Thus, T E 2

( ) = 0 in Equation (5) and a solution of the resulting equation can be given in the form 

E 

( ) r,z = A0

1+ iαd z

exp − kαd

2

⎛

⎝

⎜ ⎞

⎠

⎟ r

2

1+ iαd z

⎡

⎣

⎢ ⎤

⎦

⎥ (44) 

where A0 is the initial maximum amplitude of the wave in the output aperture and 

αd ≡ α1 + iα2 , α1 ≡ 2

kW0

2 , α2 ≡ 1

F (45) 

with W0 being the waist size of the beam as it leaves the aperture (the effective radius is aeff ( ) 0 = W0 2

) and F being the focal length of the beam. (The beam parameters αd, α1, and α2 are not to be confused 

with the absorption coefficient α of Equation (7).) Hence, one has 

E 2 = A0

2

Λ( )z

exp − kα1

Λ( )z

r

2 ⎡

⎣

⎢ ⎤

⎦

⎥ (46) 

where 

Λ( )z ≡ 1− α2 ( )z

2

+ α1

2

z

2 (47) 

and from Equation (10), 

γ 00 = π

W0

2

2

A0

2 (48) 

which is the total energy contained within the beam. Employing Equation (16) gives 
NASA/TM—2012-217634 10



S⊥ ( )z = − E 2 α2 − α1

2 + α22 ( )z

Λ( )z

⎛

⎝

⎜

⎜

⎞

⎠

⎟

⎟

( ) xxˆ + yyˆ (49) 

Equation (30) then yields 



A = 0 (50) 

and from Equations (34) and (41) 

B z( ) = − 2

kα1

⎛

⎝

⎜

⎞

⎠

⎟ α2 − α1

2 + α22 ( ) ( )z (51) 

C =

2 α1

2 + α22 ( )

kα1

 (52) 

Using these values in Equation (43) finally gives for the effective radius of the beam in a transverse plane 

at a distance z from the source aperture 

reff

2 ( )z = W0

2

2

1− α2 ( )z

2

+ α1

2

z

2 ⎡

⎣ ⎤⎦ = W0

2

2

Λ( )z (53) 

which is the value found using from Equation (46) by equating the exponential argument to -1. 

6.0 Application of the Foregoing to High Energy Laser Beam Propagation 

Through the Atmosphere 

The displacement and widening of a very high power laser beam due to atmospheric thermal 

nonlinearities are two very important performance parameters in the case of transmission of, e.g., laser 

energy from the Earth’s surface to power a spacecraft in a low Earth orbit (LEO) for subsequent 

propulsion into, e.g., a geosynchronous orbit (GEO). In this instance, a thin screen approximation can be 

made whereby only a very small length L of the entire longitudinal propagation path distance D is 

intersected by the atmosphere, i.e., L << D. To simplify matters, consider instead of the field E, its 

amplitude and phase representation 

E 

( ) r,z = e ( ) r,z exp −iφ  ⎡ ( ) r,z ⎣ ⎤⎦ (54) 

Substituting this into Equation (5) and separating the real and imaginary portions of the resulting equation 

gives 

2k ∂φ

∂z

+

∇⊥

2

e

e

+ 

( ) ∇⊥φ

2

+ k2εT T E 2

( ) = 0 (55) 

2k ∂e

∂z

− 2



∇⊥φ⋅ ∇⊥e − ∇⊥

2 ( ) φ e = 0 (56) 

Similarly, using Equation (54) in Equation (16) gives 
NASA/TM—2012-217634 11



S⊥ = − e

2 

∇φ

k (57) 

and using this intermediate result in Equations (30), (34), and (41), one has 



A = − 1

γ 00

1

k

⎛

⎝

⎜ ⎞

⎠

⎟ e

2

−∞

∞

∫ −∞

∞

∫ 

( ) ∇⊥φ dxdy (58) 

B = − 2

γ 00

1

k

⎛

⎝

⎜ ⎞

⎠

⎟ e

2 

r⊥ ⋅

 ( ) ∇⊥φ

−∞

∞

∫ −∞

∞

∫ dxdy (59) 

C = 2

γ 00

⎛

⎝

⎜

⎞

⎠

⎟

1

k

⎛

⎝

⎜ ⎞

⎠

⎟

2 

∇⊥ ( ) e

2

+ e

2 

( ) ∇⊥φ

2

− k2εTG E 2

( ) ⎡

⎣

⎢ ⎤

⎦

⎥ −∞

∞

∫ −∞

∞

∫ dxdy (60) 

where from Equation (10), 

γ 00 = e

2

−∞

∞

∫ −∞

∞

∫ dxdy (61) 

The system of Equations (55) and (56) can be considerably simplified by using the fact that the layer 

through which the laser beam will encounter thermal nonlinearities is very thin with respect to the entire 

propagation path. Hence, letting only the phase of the field be perturbed within this thin layer and, thus, 

neglecting amplitude effects, one must solve the phase equation across this thin atmospheric layer, i.e., 

2k ∂φ

∂z

+ 

( ) ∇⊥φ

2

+ k2εT T E 2

( ) = 0 (62) 

where the perturbing temperature field T E 

( ) r,z

2

( ) is given by Equation (A.3) of Appendix A. As 

shown in Appendix B, one has in this “thin screen” approximation 

φ = k x2 + y2 ( )

2Λ( )z

α2 − α1

2 + α22

( ) ( )z + πΦ( ) P0 exp −kα1y2 ( ) 1+ erf kα1 ( ) ( ) x (63) 

where 

Φ( ) P0 ≡ − P0

PT

L

2kW0

2

⎛

⎝

⎜

⎞

⎠

⎟ > 0 (64) 

which convolves the initial power P0 within the laser beam and the power PT that is characteristic of 

atmospheric thermal effects, viz, 

PT ≡ 2πρVCp

k2W0 εT α

 (65) 

The negative sign occurs in the definition of Equation (64) since, as noted earlier, εT < 0. 
NASA/TM—2012-217634 12

The evaluation of Equations (58)-(60) commences with Equation (61); by using the first relation of 

Equation (B.3) of Appendix B evaluated at the output plane of the laser in Equation (61) and performing 

the indicated integrals in plane polar coordinates gives the same result as Equation (48), i.e., 

γ 00 = πA0

2

kα1

 (66) 

The calculation of the vector 

A requires the evaluation of the transverse gradient 

∇⊥φ . Using Equation 

(63) one has 



∇⊥φ =

k α2 − α1

2 + α22 ( ) ( )z

Λ( )z

( ) xxˆ + yyˆ

+2 πΦ( ) P0 xˆ kα1

π

exp −kα1x

2 ( ) − y kˆ α1y 1+ erf kα1 ( ) ( ) x

⎡

⎣

⎢ ⎤

⎦

⎥ exp −kα1y

2 ( )

 (67) 

Substituting this expression into Equation (58) as well as that of Equation (B.3) evaluated at z = 0 (since 

A = const by Equations (30) and (31) and so can be evaluated at any longitudinal location), one finds that 

all terms involving first powers of x or y vanish upon integration in the transverse plane leaving 



A = − α1

k

Φ( ) P0 xˆ (68) 

Hence, by Equation (32), the beam is deflected by 



Δ⊥ = −z

α1

k

Φ( ) P0 xˆ (69) 

indicating that the beam moves against the atmospheric wind direction. This phenomenon is easily 

explained by the fact that the beam wants to move into regions of higher refractive index, i.e., into cooler 

regions about the beam axis. 

Since the beam deflection adds an asymmetry of the beam along the x-axis, it proves to be beneficial 

to separate out the x and y contributions of the remaining B and C parameters of Equations (59) and (60). 

To this end, Equation (59) can be written 

B = Bx + By (70) 

where 

Bx = − 2

γ 00

1

k

⎛

⎝

⎜ ⎞

⎠

⎟ e

2

xxˆ ⋅

 ( ) ∇⊥φ

−∞

∞

∫ −∞

∞

∫ dxdy (71) 

By = − 2

γ 00

1

k

⎛

⎝

⎜ ⎞

⎠

⎟ e

2

yyˆ⋅

 ( ) ∇⊥φ

−∞

∞

∫ −∞

∞

∫ dxdy (72) 

Using Equation (67) and once again noting that terms containing first powers of x and y vanish upon 

integration, one obtains 
NASA/TM—2012-217634 13

Bx = − α2

kα1

, By = − α2

kα1

+

1

k

π

2

Φ( ) P0 (73) 

One must now separate Equation (60) into x and y component contributions. Thus, without loss of 

generality, one can write 

C = Cx + Cy (74) 

where 

Cx = 2

γ 00

⎛

⎝

⎜

⎞

⎠

⎟

1

k

⎛

⎝

⎜ ⎞

⎠

⎟

2

xˆ ⋅



∇⊥ ( ) e

2

+ e

2

xˆ ⋅

 ( ) ∇⊥φ

2 ⎡

⎣

⎢ ⎤

⎦

⎥ −∞

∞

∫ −∞

∞

∫ dxdy

⎧

⎨

⎪

⎩⎪

− Q k2εTG E 2

( ) ⎡

⎣ ⎤⎦ −∞

∞

∫ −∞

∞

∫ dxdy

⎫

⎬

⎪

⎭⎪

 (75) 

and 

Cy = 2

γ 00

⎛

⎝

⎜

⎞

⎠

⎟

1

k

⎛

⎝

⎜ ⎞

⎠

⎟

2

yˆ⋅



∇⊥ ( ) e

2

+ e

2

yˆ⋅

 ( ) ∇⊥φ

2 ⎡

⎣

⎢ ⎤

⎦

⎥ −∞

∞

∫ −∞

∞

∫ dxdy

⎧

⎨

⎪

⎩⎪

− ( ) 1− Q k2

εTG E 2

( ) ⎡

⎣ ⎤⎦ −∞

∞

∫ −∞

∞

∫ dxdy

⎫

⎬

⎪

⎭⎪

 (76) 

where the G E 2

( ) terms have been separated out for convenience. It is very important to note that 

G E 2

( ) cannot be arbitrarily apportioned between the arbitrary x and y components of the coefficient C. 

The G E 2

( ) term has its origins with the conservation of energy given by Equation (15) and is thus 

associated with the deposition of energy within the propagating laser beam. The beam deflection adds 

asymmetry to the beam and energy will be deposited along the x and y directions differently. In what is to 

follow, a fraction Q of the energy density will be taken to correspond to that along the x axis and the 

remaining 1 – Q fraction along the y axis. An approximation procedure will be given below to determine 

the division of energy among these components and to find a value for Q. 

Consider, at the outset, the evaluation of the first integrand terms of Equations (75) and (76). Once 

again, working in the z = 0 plane, one has from the first of Equation (B.3) 

xˆ ⋅



∇⊥ ( ) e

2

= A0

2

k2α1

2

x

2 exp −kα1 x2 + y2 ( ) ( ) (77) 

and a similar relation for the y component. One thus has 

xˆ ⋅



∇⊥ ( ) e

2

−∞

∞

∫ −∞

∞

∫ dxdy = yˆ⋅



∇⊥ ( ) e

2

−∞

∞

∫ −∞

∞

∫ dxdy = A0

2

π

2 (78) 

Using Equation (67), one also has 

e

2

xˆ ⋅

 ( ) ∇φ 2

= A0

2 exp −kα1 x

2 + y2 ( ) ( ) kα2 ( ) x

2

+ 4kα2 xΦ( ) P0 kα1 ⎡

⎣

⋅exp −kα1x

2 ( ) exp −kα1y2 ( ) + 4Φ2 ( ) P0 kα1 ⋅⋅exp −2kα1x2 ( ) exp −2kα1y2 ( )⎤

⎦

 (79) 

and similarly, 
NASA/TM—2012-217634 14

e

2

yˆ⋅

 ( ) ∇φ 2

= A0

2 exp −kα1 x

2 + y2 ( ) ( ) y2 ( ) kα2

2 − 4k2

α1α2 πΦ( ) P0 ⎡

⎣

⋅exp −kα1y

2 ( ) 1+ erf kα1 ( ) ( ) x + +4k2

α1

2

πΦ2 ( ) P0 exp −2kα1y

2 ( ) 1+ erf kα1 ( ) ( ) x

2 ⎤

⎦

⎥

 (80) 

These yield 

e

2

−∞

∞

∫ −∞

∞

∫ xˆ ⋅

 ( ) ∇φ 2

dxdy = A0

2 π

2

α2

2

α1

2 + A0

2 4

3

πΦ2 ( ) P0 (81) 

and 

e

2

−∞

∞

∫ −∞

∞

∫ yˆ⋅

 ( ) ∇φ 2

dxdy = A0

2 π

2

α2

2

α1

2 − A0

2

π

π

2

α2

α1

Φ( ) P0 + A0

2

π3/2 8

9

π

3

Φ2 ( ) P0 (82) 

The evaluation of the second integrand term in Equation (76) commences with the use of 

Equation (40), i.e., 

k2εT G E 2

( ) −∞

∞

∫ dx dy =

−∞

∞

∫ 1

2

k2εT T E 2

−∞

∞

∫ dx dy

−∞

∞

∫ (83) 

One now employs Equations (46) and (A.3) in the right side of this relation and obtains in the z = 0 plane 

T E 2

−∞

∞

∫ dxdy z=0 =

−∞

∞

∫ A0

4

2 β π

kα1

exp −kα1x

2 ( ) −∞

∞

∫ −∞

∞

∫ exp −2kα1y

2 ( ) 1+ erf kα1 ( ) ( ) x dxdy

= A0

4

2 β 1

kα1

⎛

⎝

⎜

⎞

⎠

⎟

3/2

π3/2

2

 (84) 

Using this in Equation (83) and rearranging terms gives 

k2εT G E 2

( ) −∞

∞

∫ dx dy =

−∞

∞

∫ P0

PT

πA0

2

4

⎛

⎝

⎜

⎞

⎠

⎟

π

2 (85) 

It may be noted that Equation (85) is independent of the thickness L of the absorbing medium. This is due 

to the level of estimation that the ‘thin screen’ approximation yields; in that approximation, only the 

phase and not the amplitude of the field intensity of the beam is perturbed. It would only be with the next 

approximation beyond that used in Appendix B that would add a perturbation to the amplitude of the 

field. This limitation must be kept in mind in the applications. 

As discussed above, the division of this quantity among the x and y components cannot be arbitrarily 

done. The value of Q, introduced above, will be calculated shortly. Thus, substituting Equations (66), 

(78), (81) and Equation (85) into Equation (75) results in 

Cx = α1

k

1+

α2

2

α1

2 +

8

3

Φ2 ( ) P0

⎛

⎝

⎜

⎞

⎠

⎟ + Q π

2

2

kL

⎛

⎝

⎜ ⎞

⎠

⎟ Φ( ) P0 (86) 

and, similarly, putting Equations (66), (78), (82), and Equation (85) into Equation (76) gives 
NASA/TM—2012-217634 15

Cy = α1

k

1+

α2

2

α1

2 − 2 α2

α1

π

2

Φ( ) P0 +

2

3

π

3

Φ2 ( ) P0

⎛

⎝

⎜

⎞

⎠

⎟ + ( ) 1− Q π

2

2

kL

⎛

⎝

⎜ ⎞

⎠

⎟ Φ( ) P0 (87) 

In applying Equation (43) to calculate the beam radius about the local axis, care must be taken in the 

specifications for the B and C parameters using Equations (73), (86) and (87). In fact, the derivation of 

Equations (73), (86) and (87) require the corresponding parameters in Equation (43) to be multiplied by a 

factor of 2, i.e., 

reff , x

2 ( )z = 1

2

( ) 2Cx z

2 + ( ) 2Bx ( ) 0 z + aeff

2 ( ) 0 − 

( ) A⋅ xzˆ 2

 (88) 

and 

reff , y

2 ( )z = 1

2 ( ) 2Cy z

2 + ( ) 2By ( ) 0 z + aeff

2 ( ) 0 − 

( ) A⋅ yzˆ 2

 (89) 

Using Equations (88) and (89), one has from the foregoing, after rearrangement of factors and 

simplification, 

α1kreff , x

2 ( )z = reff , x

2

W0

2 = 1− α2 ( )z

2

+ α1

2 1+

5

3

Φ2 ( ) P0

⎛

⎝

⎜ ⎞

⎠

⎟ z

2 + Qα1

π

2

2

L

⎛

⎝

⎜ ⎞

⎠

⎟ Φ( ) P0 z

2 (90) 

and 

α1kreff , y

2 ( )z = reff , y

2

W0

2 = 1− α2 ( )z

2

+ α1

2 1+

16

9

π

3

Φ2 ( ) P0

⎛

⎝

⎜ ⎞

⎠

⎟ z

2

+( ) 1− Q α1

π

2

2

L

⎛

⎝

⎜ ⎞

⎠

⎟ Φ( ) P0 z

2 + 2α1

π

2

Φ( ) P0 z 1− α2 ( )z

 (91) 

The case of a focused beam wave, i.e., F = D will be considered, without addressing the methodology 

with which to achieve it in practice. In this instance, the above relations reduce to 

reff , x

2

W0

2

foc

= α1

2 1+

5

3

Φ2 ( ) P0

⎛

⎝

⎜ ⎞

⎠

⎟ D2 + Qα1

π

2

2

L

⎛

⎝

⎜ ⎞

⎠

⎟ Φ( ) P0 D2 (92) 

and 

reff , y

2

W0

2

foc

= α1

2 1+

16

9

π

3

Φ2 ( ) P0

⎛

⎝

⎜ ⎞

⎠

⎟ D2 + ( ) 1− Q α1

π

2

2

L

⎛

⎝

⎜ ⎞

⎠

⎟ Φ( ) P0 D2 (93) 

An expression for the value of Q can now be obtained. The first member of the right side of Equation (92) 

is the semi-axis of the beam in the x direction and that of Equation (93) is the semi-axis of the beam in the 

y direction. As a first approximation, it can be assumed that the ratio of the lengths of these semi-axes is 

equal to the corresponding ratio of the energy deposited along these axes, i.e., 
NASA/TM—2012-217634 16

1+

5

3

Φ2 ( ) P0

1+

16

9

π

3

Φ2 ( ) P0

= Q

1− Q (94) 

Solving this expression for Q yields 

Q = Q P( ) 0 =

PT

2 + K1P02 Lα1

4

2PT

2 + ( ) K1 + K2 P02 Lα1

4

 (95) 

where K1 = 5 3 and K2 = 16π 9 3 . This expression is well behaved at the limits lim

P0 →0

Q = 1 2

(indicating that the energy is deposited equally along the axes when no thermal action is prevalent) and 

lim

P0 →∞Q = 0.341 . Once again, it must be remembered that this approximation procedure is in lieu of a 

second-order solution to the equations of propagation in which the field amplitude of the beam wave is 

also given as a function of propagation conditions along the beam. 

Numerical application of the foregoing to focused beam propagation in a nonlinear thermal medium 

commences with the rewriting of Equations (92) and (93) in a dimensionless form. Additionally, the 

expression for the corresponding beam deflection 

Δ⊥ as given by Equation (69) will be treated. First, 

dealing with the latter expression, Equation (69) and be rewritten in the dimensionless form (deflection 

from intended axis normalized by the initial beam radius) 

Δn ≡



Δ⊥

W0

= 2 D

kW0

2 Φ( ) P0 = 1

2

ldΦ( ) P0 (96) 

where ldD ≡ 2D kW0

2 is the diffraction length associated with the aperture size and the total propagation 

distance D to the target. Figures 1 and 2 depict the variation of the normalized deflection of the beam 

versus the beam power in kilowatts, for ldD = 0.02 and ldD = 0.2 , respectively; these correspond to 

realistic beam and propagation conditions. The direction of the deflection is along the –x axis, i.e., it is 

directed into the wind where regions of higher refractive index prevail. As heat is removed from the 

propagation channel along the x axis, it continues to accumulate along the y axis of the channel. This 

gives rise to an asymmetry in the beam cross-section as is shown by considering the beam radii. 

Similarly, for the beam radii at the intended focus, one can write Equation (92) as 

reff , x

2 W02

( ) α1D 2

foc

≡ rF, x

2 = 1+

5

3

Φ2 ( ) P0

⎛

⎝

⎜ ⎞

⎠

⎟ + Q π

2

2

α1L

⎛

⎝

⎜

⎞

⎠

⎟ Φ( ) P0 (97) 

The left side is essentially the square of the ratio of the beam radius in the x direction normalized to the 

intended radius of the focal spot on the target, 2D kW0 . The same for Equation (93), i.e., 

rF, y

2 = 1+

16

9

π

3

Φ2 ( ) P0

⎛

⎝

⎜ ⎞

⎠

⎟ + ( ) 1− Q π

2

2

α1L

⎛

⎝

⎜

⎞

⎠

⎟ Φ( ) P0 (98) 
NASA/TM—2012-217634 17

In both of these relations, the quantity ldL ≡ α1L is the diffraction length with respect to the aperture size 

and the thickness of the atmospheric layer L. Figures 3 and 4 depict the beam radii variation versus beam 

power for two representative values of ldL . 
NASA/TM—2012-217634 18

In addition to the evaluation of the deleterious effects of the atmosphere on the propagation of highenergy laser beams, the methodology introduced above is amenable to the analysis of the design and 

specification of adaptive optics approaches whereby the phase front of the transmitted laser beam is 

appropriately modified to mitigate the atmospheric effects (Ref. 10). For example, due to the elliptical 

profile of the radiation in the transverse plane, the use of the well-known Zernike polynomials (Ref. 11), 

which are defined on a unit circle, are not appropriate for this application. Instead, the transmitted phase 

front φT ( ) x, y at the output aperture of the transmitter can be expanded in a set of orthogonal 
NASA/TM—2012-217634 19

polynomials associated with a Gaussian weight function, i.e., the Hermite polynomials Hn(x) and Hm(y). 

Thus, 

φT ( ) x, y = anm

m

N

∑

n

M

∑ Hn ( ) x Hm ( ) y (99) 

The expansion coefficients anm are determined by applying one of several ‘performance metrics’. For 

example, for power beaming applications, it is desired to shape the phase front so as to minimize the 

radius of the beam at the target. Such a representation is straightforward to use in Equation (63) and 

subsequently, in Equations (71), (72), (75) and (76) to obtain the corrected radius values for the beam. 

7.0 Summary and Conclusion 

The structure of the nonlinear parabolic equation that describes the propagation of a high energy laser 

beam through a medium with thermal nonlinearities is exploited to describe the moments of the 

associated electric field in order to capture the universal behavior of the displacement and radius of the 

beam; these are two of the more significant performance parameters for the operation of the high energy 

laser transmission system. Once this was done, a simple thin phase screen propagation model of the Earth 

atmosphere with a wind was advanced and applied to the equations for the beam displacement and radius. 

The presence of the wind of course broke the symmetry of the beam into one with an elliptical cross 

section. Expressions were then developed for the behavior of the semi-axes of the beam for the 

propagation path out of the atmosphere. Graphical results were then presented for the deleterious action of 

the atmosphere on the propagation of the beam. Since only phase perturbations were admitted into the 

model, all the effects placed on the propagating beam are phase based and can thus be mitigated at the 

transmitter. The elliptical cross-section of the beam dictates that the circular Zernike polynomials are not 

sufficient to describe the phase front and that a Hermite polynomial based phase decomposition be used. 

Such a model was exploited for the analysis of laser power beaming for space propulsion (Ref. 12). 

NASA/TM—2012-217634 21

Appendix A 

Equation (7) can be simplified straight away by first specializing the wind velocity to be along the x

axis of the coordinate system, 

V = Vxˆ , where xˆ is the unit vector. Further, the diffusivity ∇2

T will be 

taken to have contributions only in the direction transverse to the laser beam propagation, i.e., 

∇2T ≈ ∇⊥

2 T . Finally, only the stationary heating case will be considered for this particular discussion 

whereby the time derivative can be dropped. (That is, the CW radiation is taken to be acting long enough 

for the steady state case to be achieved.) Applying these considerations to Equation (7) allows it to be 

written 

V ∂T

∂x

− χ∇⊥

2 T = αc

8πρCp

E 2 (A.1) 

where χ≡κ ρCp is the associated thermal diffusivity. There are two heat transfer mechanisms that are 

described by this relation, i.e., forced convection due to the wind given by the first term on the left side of 

the equation and diffusion by the second term on the left side. The order of magnitude of the 

corresponding derivatives are ∂T ∂x ~ T reff and ∇⊥

2 T ~ T reff2 where reff is the effective laser beam 

radius. Consider now the ratio of the coefficients of the two terms within the brackets of Equation (A.1), 

i.e., η ≡ Vreff χ . This dimensionless ratio is known in fluid mechanics as the Peclet number. Taking the 

nominal wind velocity V = 4.4 m/s (10 mi/hr), reff = 10 m and using the documented value for the thermal 

diffusivity of the atmosphere χ = 2.12 ×10−5 m2

/s, one has that η = 2.1×106 >> 1 allowing one to 

neglect the second term in the brackets. In this instance, Equation (A.1) becomes 

∂T

∂x = β E 2 (A.2) 

where β≡αc ( ) 8πVρCp . This equation can now be easily integrated using Equation (46). The result is 

T E 

( ) r,z

2

( ) = β E( ) x′, y,z

2

−∞

x

∫ dx′

= β A0

2

2

π

kα1Λ( )z

⎛

⎝

⎜

⎞

⎠

⎟

1/2

exp − kα1

Λ( )z

y

2 ⎡

⎣

⎢ ⎤

⎦

⎥ 1+ erf kα1

Λ( )z

⎛

⎝

⎜

⎞

⎠

⎟

1/2

x

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

⎛

⎝

⎜

⎜

⎞

⎠

⎟

⎟

 (A.3) 

Now in order to substantiate Equation (40) for this case, one proceeds to define 

IT Ex( ) 2

( ),G Ex( ) 2

( ) ( ) ≡ T E 2 − 2G E 2

( ) ⎡

⎣ ⎤⎦ −∞

∞

∫ dx dy

−∞

∞

∫ (A.4) 

and attempts to show that the expression is equal to zero. To this end, one forms the derivative 

dI

dx = ∂T

∂x

E 2+

∂ E 2

∂x

T − 2 ∂G

∂x

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ −∞

∞

∫ dx dy

−∞

∞

∫ (A.5) 
NASA/TM—2012-217634 22

which, if also zero, can be used to argue that Equation (A.4) vanishes do to conditions at x → ∞ . From 

Equation (13), 

G( ) Ξ( ) x = T ( ) Ξ′

0

Ξ( ) x

∫ dΞ′, Ξ( ) x ≡ E x( ) 2

 (A.6) 

Differentiating Equation (A.6), 

∂G

∂x = ∂Ξ( ) x

∂x

T ( ) Ξ +

∂T ( ) Ξ′

∂x

0

Ξ

∫ dΞ′ (A.7) 

substituting into Equation (A.5), using Equation (A.2), performing the indicated functional integration, 

and simplifying finally gives 

dI

dx = − ∂ E 2

∂x

T

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ −∞

∞

∫ dx dy

−∞

∞

∫ (A.8) 

which does not appear to vanish. To be sure, proceeding to integrate Equation (A.8) by parts yields 

dI

dx = E 2 ∂T∂x

⎡

⎣

⎢ ⎤

⎦

⎥ −∞

∞

∫ dx dy

−∞

∞

∫ = −β E 4

−∞

∞

∫ dx dy

−∞

∞

∫ (A.9) 

upon remembering Equation (A.2). The E 4 dependency is reminiscent of the Kerr effect that appears at 

field intensities much higher than for thermal blooming effect. To be sure, in order to secure numerical 

estimates, one uses Equation (46) in Equation (A.9) to obtain 

dI

dx = βA0

4 π

4

W0

2 (A.10) 

To rid of the amplitude factors and put them into a more familiar quantity, i.e., the total power P0 within 

the laser beam, one uses the relation 

P0 ≡ cA0

2

8π

⎛

⎝

⎜

⎞

⎠

⎟

πW0

2

2

⎛

⎝

⎜

⎞

⎠

⎟ (A.11) 

as well as remember the definition of β to write Equation (A.10) as 

dI

dx = 16αP0

2

cW0

2

VρCp

 (A.12) 

Using the values of the parameters given above and additionally taking W0 ≈ 1 , this gives 

dI dx ≈ 4 ×10−5 for P0 = 1×106 watts. Hence, dI dx effectively vanishes for nominal laser power and 

thus Equation (A.4) can also be taken to vanish, i.e., 
NASA/TM—2012-217634 23

T E 2 − 2G E 2

( ) ⎡

⎣ ⎤⎦ −∞

∞

∫ dx dy

−∞

∞

∫ ≈ 0 (A.13) 

A more complete treatment would involve a more general solution to the temperature diffusion equation, 

Equation (A.1), but will not be dealt with here. 

NASA/TM—2012-217634 25

Appendix B 

Writing the total phase as the sum of the unperturbed phase φ0 (i.e., defined by Equation (62) with 

T E 2

( ) = 0 ) and the perturbed component φ , one has, upon introducing φ = φ0 + φ into Equation (62), 

the system of equations 

2k ∂φ0

∂z

+ 

( ) ∇⊥φ0

2

= 0 (B.1) 

2k ∂φ

∂z

+ 2



∇⊥φ0 ⋅



∇⊥φ + 

∇⊥φ ( )2

+ k2εT T E 2

( ) = 0 (B.2) 

Now, equating Equations (44) and (54), one finds that 

e 

( ) r,z = A0

1+ iαd z

exp − kα1

2Λ( )z

r

2 ⎛

⎝

⎜

⎞

⎠

⎟

φ0



( ) r,z = kr2

2Λ( )z

α2 − α1

2 + α22 ( ) ( )z

 (B.3) 

the second relation of which, indeed, can be shown to satisfy Equation (B.1). 

A solution for Equation (B.2) commences with multiplying it by 1/ k and considering k → ∞

(relative to all the other inverse lengths in the problem) which gives 

2 ∂φ

∂z

+ kεT T E 2

( ) = 0

 (B.4) 

One needs to use Equation (A.3) to integrate this expression. To this end, one applies the following ‘thin 

screen’ approximation; For the case considered here, the perturbing atmospheric layer of thickness L, in 

which the thermal nonlinearities take place, is such that L << D where D is the total distance to the target 

between LEO and GEO. Thus, within this region, relative to the total propagation distances involved, one 

can let Λ( )z ≈ 1 for the longitudinal distances z occurring within the atmospheric layer. Thus, one can 

write 

T E 2

( ) ≈ βA0

2 π

W0

2

exp −kα1y

2 ⎡

⎣ ⎤⎦

1

2

⎛

⎝

⎜ ⎞

⎠

⎟ 1+ erf kα1 ⎡ x ⎣ ⎤ ( )⎦ (B.5) 

Substituting this result into Equation (B.4) and integrating across the thickness L of the atmosphere 

yields, upon rearranging factors, 

φ = − P0

PT

⎛

⎝

⎜

⎞

⎠

⎟

πL

kW0

2

⎛

⎝

⎜

⎞

⎠

⎟

1

2

⎛

⎝

⎜ ⎞

⎠

⎟ exp −kα1y

2 ( ) 1+ erf kα1 ( ) ( ) x (B.6) 

where 

P0 ≡ cA0

2

8π

⎛

⎝

⎜

⎞

⎠

⎟

πW0

2

2

⎛

⎝

⎜

⎞

⎠

⎟ (B.7) 
NASA/TM—2012-217634 26

is the initial power of the beam wave of waist radius W0 and 

PT ≡ 2πρVCp

k2W0 εT α

 (B.8) 

is the characteristic power for the onset of thermal blooming. Defining 

Φ( ) P0 ≡ − P0

PT

L

2kW0

2

⎛

⎝

⎜

⎞

⎠

⎟ > 0 (B.9) 

(since εT < 0 ), Equations (B.3) and (B.6) yield for the total phase of the laser beam leaving the thin 

atmospheric layer with a wind of velocity V

φ =

k x2 + y

2 ( )

2Λ( )z

α2 − α1

2 + α22 ( ) ( )z + πΦ( ) P0 exp −kα1y

2 ( ) 1+ erf kα1 ( ) ( ) x (B.10) 
NASA/TM—2012-217634 27

References 

1. R. M. Manning, Stochastic Electromagnetic Image Propagation Through the Atmosphere (McGrawHill Book Co., New York, NY, 1993), Chapter 2. 

2. For a good recent treatment, see A. D. Wheelon, Electromagnetic Scintillation: Volume 2, Weak 

Scattering (Cambridge University Press, 2003), and references therein. 

3. S. A. Akhmanov, D. P. Krindach, A. V. Migulin, A. P. Sukhorukov, and R. V. Khokhlov, “Thermal 

Self-Action of Laser Beams,” IEEE J. Quantum Electron. QE-4 (10), pp. 568-575 (1968). 

4. D. C. Smith, “High-Power Laser Propagation: Thermal Blooming,” Proc IEEE 65 (12), pp. 1679-

1714 (1977). 

5. V. V. Vorob’ev, “Thermal Blooming of Laser Beams in the Atmosphere,” Prog. Quantum Electron.,

15, pp. 1-152 (1991). 

6. M. V. Goldman and D. R. Nicholson, “Virial Theory of Direct Langmuir Collapse,” Phys. Rev. Lett. 

41 (6), pp. 406-410 (1978). 

7. D. W. McLaughlin, G. C. Papanicolaou, C. Sulem, and P. L. Sulem, “Focusing Singularity of the 

Cubic Schrodinger Equation,” Phys. Rev. A 34 (2), pp. 1200-1210 (1986). 

8. A. Zh. Muradyan, “Resonance Self-Focusing in the Limit of Strong Nonlinearity Saturation,” Sov. 

Phys. JETP 65 (6), pp. 1111-1114 (1987). 

9. M. V. Goldman, K. Rypdal, and B. Hafizi, “Dimensionality and Dissipation in Langmuir Collapse,” 

Phys. Fluids 23 (5), pp. 945-955 (1980). 

10. L. C. Bradley and J. Herrmann, “Phase Compensation for Thermal Blooming,” Appl. Opt. 13 (2), pp. 

331-334 (1974). 

11. D. E. Novoseller, “Zernike-Ordered Adaptive-Optics Correction of Thermal Blooming,” J. Opt. Soc. 

Am. A5 (11), pp. 1937-1942 (1988). 

12. P. George and R. Beach, “Beamed-Energy Propulsion (BEP) Study,” NASA Technical 

Memorandum, NASA/TM—2012-217014, February 2012, Appendix G. 
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-07-2012 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

High Energy Laser Beam Propagation in the Atmosphere: The Integral Invariants of the 

Nonlinear Parabolic Equation and the Method of Moments 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Manning, Robert, M. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 439432.04.07.05 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-18279 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2012-217634 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 36, 44, 70, and 74 

Available electronically at http://www.sti.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

The method of moments is used to define and derive expressions for laser beam deflection and beam radius broadening for high-energy 

propagation through the Earth’s atmosphere. These expressions are augmented with the integral invariants of the corresponding nonlinear 

parabolic equation that describes the electric field of high-energy laser beam to propagation to yield universal equations for the 

aforementioned quantities; the beam deflection is a linear function of the propagation distance whereas the beam broadening is a quadratic 

function of distance. The coefficients of these expressions are then derived from a thin screen approximation solution of the nonlinear 

parabolic equation to give corresponding analytical expressions for a target located outside the Earth’s atmospheric layer. These equations, 

which are graphically presented for a host of propagation scenarios, as well as the thin screen model, are easily amenable to the phase 

expansions of the wave front for the specification and design of adaptive optics algorithms to correct for the inherent phase aberrations. This 

work finds application in, for example, the analysis of beamed energy propulsion for space-based vehicles. 

15. SUBJECT TERMS

Electromagnetic wave transmission; Wave propagation; High power lasers; Laser beams 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

34 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18