Spectral Solution of the Helmholtz and Paraxial Wave 

Equations and Classical Diffraction Formulae 

by Timothy M. Pritchett 

ARL-TR-3179 March 2004 

Approved for public release; distribution unlimited. 
NOTICES 

Disclaimers 

The findings in this report are not to be construed as an official Department of the Army position unless 

so designated by other authorized documents. 

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the 

use thereof. 

Destroy this report when it is no longer needed. Do not return it to the originator. 
Army Research Laboratory 

Adelphi, MD 20783-1197 

ARL-TR-3179 March 2004 

Spectral Solution of the Helmholtz and Paraxial Wave 

Equations and Classical Diffraction Formulae 

Timothy M. Pritchett 

Sensors and Electron Devices Directorate, ARL 

Approved for public release; distribution unlimited. 
ii

REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188

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 information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the 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)

March 2004 

2. REPORT TYPE 

Final 

3. DATES COVERED (From - To) 

November 2003 to December 2003 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

4. TITLE AND SUBTITLE 

Spectral Solution of the Helmholtz and Paraxial Wave Equations and Classical Diffraction 

Formulae 

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER 

5e. TASK NUMBER 

6. AUTHOR(S) 

Timothy M. Pritchett 

5f. WORK UNIT NUMBER 

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

U.S. Army Research Laboratory 

Attn: AMSRD-ARL-SE-EO 

2800 Powder Mill Road 

Adelphi, MD 20783-1197

8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

ARL-TR-3179 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S) 

U.S. Army Research Laboratory 

2800 Powder Mill Road 

Adelphi, MD 20783-1197 11. SPONSOR/MONITOR'S REPORT 

 NUMBER(S) 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Approved for public release; distribution unlimited. 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

This report examines the general nonlinear vector wave equation implied by the Maxwell equations in a nonmagnetic, isotropic 

medium and discusses the various approximations under which this general result reduces to the familiar scalar Helmholtz 

equation and the paraxial wave equation. We see why the Helmholtz equation may be regarded as a singular perturbation of the 

paraxial wave equation and how some of the difficulties arising in the solution of the former partial differential equation are

related to this fact. Standard integral transform methods are used to obtain general solutions of the Helmholtz equation in a 

linear medium and of the paraxial wave equation in a linear medium. We show that these solutions are equivalent, respectively,

to the exact Rayleigh-Sommerfeld diffraction integral and the Rayleigh-Sommerfeld integral in the Fresnel approximation. We 

discuss the use of these linear solutions in numerical procedures for treating the corresponding nonlinear beam propagation 

problems. The general linear solutions are specialized to situations with axial symmetry, and the results are used to treat the 

example of a clipped Gaussian beam. 

15. SUBJECT TERMS

Diffraction, Fresnel approximation, Helmholtz equation 

16. SECURITY CLASSIFICATION OF: 19a. NAME OF RESPONSIBLE PERSON

Timothy M. Pritchett 

a. REPORT 

UNCLASSIFIED 

b. ABSTRACT 

UNCLASSIFIED

c. THIS PAGE 

UNCLASSIFIED

17. LIMITATION 

OF ABSTRACT 

UL 

18. NUMBER 

OF PAGES 

28 19b. TELEPHONE NUMBER (Include area code) 

(301) 394-0332 

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18 
iii

Contents 

List of Figures iv

1. Introduction 1

2. Scalar Electromagnetic Wave Equations and Their Validity 2

2.1 General Wave Equation for Nonlinear Optical Media....................................................2 

2.2 Homogeneous Helmholtz Equation.................................................................................4 

2.3 Paraxial Wave Equation ..................................................................................................4 

3. A Précis of Classical Scalar Diffraction Theory 5

3.1 History.............................................................................................................................5 

3.2 Rayleigh-Sommerfeld Diffraction Integral .....................................................................5 

3.3 Preliminary Approximations ...........................................................................................7 

3.4 Fresnel Approximation....................................................................................................7 

3.5 Fraunhofer Limit .............................................................................................................7 

4. Solution of the Paraxial Wave Equation 8

4.1 Helmholtz Equation as a Singular Perturbation of the Paraxial Wave Equation ............8 

4.2 Convolution Kernel for the Paraxial Wave Equation......................................................9 

4.3 A Remark on Convergence ...........................................................................................10 

4.4 Numerical Computation of the Fresnel Convolution Integral.......................................10 

5. Solution of the Helmholtz Equation 11

6. Axial Symmetry 12

6.1 Cylindrically Symmetric Restriction of the General Fresnel Kernel ............................12 

6.2 Solution of the Paraxial Wave Equation Via Hankel Transform ..................................13 

6.3 Cylindrically Symmetric Fraunhofer Kernel.................................................................14 

7. Example: Clipped Gaussian Beam 14

7.1 Fresnel Diffraction of a Clipped Gaussian Beam..........................................................15 

7.2 Fraunhofer Diffraction of a Clipped Gaussian Beam....................................................15 
iv

7.3 Limiting Case: Top Hat Beam (Uniformly Illuminated Aperture) ..............................15 

8. Summary 16

9. References 18

Distribution 21 

List of Figures 

Figure 1. Diffraction by an aperture: geometry of the problem......................................................6 
1

1. Introduction 

To accurately describe the propagation of a light beam through a medium, one must account, 

roughly speaking, for three phenomena: diffraction, refraction, and absorption. Diffraction 

refers to the tendency of wave phenomena such as light or sound to deviate from rectilinear 

propagation when a portion of the wave front is obstructed in some way. Diffractive effects are 

present to some degree in every real optical system, simply because no real system can be 

infinite in extent; every real system must necessarily contain an effective aperture of some sort. 

Refraction is related to the speed with which a wave travels through a medium. It, too, can result 

in the deviation of a beam from straight-line propagation, as when the refractive properties of the 

medium depend on position in such a way that different portions of the beam experience 

different amounts of refraction. Absorption is the process responsible for the change in 

amplitude of the wave as it propagates through the medium and energy is transferred from the 

wave to the medium or from the medium to the wave (gain). Both refraction and absorption 

depend on the amplitude of the propagating wave. In certain media, this dependence is highly 

nonlinear, while in other materials, such as “slow” (fluence dependent) excited state absorbers 

(1), it is, to be sure, linear in the wave amplitude but reflects the entire history of the interaction 

between the medium and the propagating beam. In general, calculating the effects of 

propagation of a high-intensity laser beam through all but the simplest media presents a very 

difficult problem. 

In these “simple” media (isotropic media exhibiting refractive and absorptive responses that are 

both linear and time independent), beam propagation essentially reduces to pure diffraction. 

However, even a pure diffraction problem in an uncomplicated geometry can generally not be 

solved without recourse to various approximations, and these impose limits on the resulting 

solutions. Nevertheless, solutions do exist. 

These solutions are very valuable. The special case of beam propagation (diffraction) in an 

isotropic, linear, and time-independent medium is extremely important in practice. What is 

more, the solution of the linear diffraction problem can be employed to attack the more general 

problem of beam propagation in a nonlinear medium through the use of a numerical procedure 

that separately computes the effects of diffraction, nonlinear refraction, and nonlinear absorption 

for each propagation step. This is the “split step” Fourier method introduced by Fleck, Morris, 

and Feit (2) (sometimes referred to simply as the Beam Propagation Method). 

The remainder of this report is organized as follows. Section 2 opens with an examination of the 

general nonlinear optical wave equation and the conditions under which it reduces to the scalar 

Helmholtz equation. This is followed by a careful derivation of the paraxial wave equation. 

Section 3 briefly reviews the salient results of classical scalar diffraction theory. In sections 4 

and 5, standard integral transform methods are used to obtain a general solution to the paraxial 
2

wave equation in a linear medium and the Helmholtz equation in a linear medium, respectively. 

These results underscore the fact (perhaps not widely appreciated) that the exact RayleighSommerfeld integral is equivalent to the scalar Helmholtz equation, while the RayleighSommerfeld integral in the Fresnel approximation is equivalent to the paraxial wave equation. 

The general results of section 4 are specialized to situations with cylindrical symmetry in section 

6. In section 7, these cylindrically symmetric results are applied to a practical example—a 

Gaussian beam clipped by an aperture. 

2. Scalar Electromagnetic Wave Equations and Their Validity 

2.1 General Wave Equation for Nonlinear Optical Media 

The derivation of an electromagnetic wave equation from the Maxwell equations is presented in 

virtually any textbook on optics or electromagnetic theory. The most general form (3) of the 

wave equation in a nonmagnetic medium containing neither free charges nor free currents is, in 

Gaussian units, 

P( , , ) 4 ( , , ) 1 ( , , ) 2

2

2 2

2

2 x z t

c t

x z t

c t

x z t r r r r r r r r

∂

∂ Ε = − ∂∂ ∇ ×∇ × Ε +π . (1) 

(In anticipation of future developments, we place one of the coordinates, z, on a separate footing 

and combine the remaining two coordinates into the two-component vector x

r

, which designates 

positions in the transverse plane.) In equation 1, P(x,z,t) r r is the total material polarization. It is 

usually expressed as a power series in the electric field ) (x,z,t r r Ε , with its ith component given 

by 

K

r r r r

r r r

r r r

′′′ ′′ ′ − ′ − ′′ − ′′′ Ε ′ Ε ′′ Ε ′′′ +

′′ ′ − ′ − ′′ Ε ′ Ε ′′ +

= ′ − ′ Ε ′ +

∫ ∫ ∫

∫ ∫

∫

∞

−∞

∞

−∞

∞

−∞

∞

−∞

∞

−∞

∞

−∞

( , , , , ) ( , , ) ( , , ) ( , , )

( , , , ) ( , , ) ( , , )

P ( , , ) ( , , ) ( , , )

(3)

(2)

(1)

dt dt dt x z t t t t t t x z t x z t x z t

dt dt x z t t t t x z t x z t

x z t dt x z t t x z t

i j k l j k l

i j k j k

i i j j

χ

χ

χ

 

In this expression, the quantities ( , , , ) ( ) K r K x z t n χ i j are components of the nth-order time-dependent 

electric susceptibility tensor, and the repeated indices j, k, l, .... carry an implied sum (Einstein 

summation convention). The convolutions in each term of the above expansion for ) P (x,z,t i

r

become simple products upon Fourier transformation. 
3

We consider now the single-frequency component of the electric field oscillating at angular 

frequency ω: 

(x,z,t) E(x,z)e c.c. i t Ε = + r − ω r r r (2) 

We employ the well-known vector identity for the curl of a curl and the constitutive relation in 

the form 

P n E r r 4 ( 1) 2 π = −

to recast equation 1 as 

( ) ( ) n x z E t n E

c

E

c

n E E r r r r r r r r 2

0

2

2

2

2

2

0

2

2 ∇ ∇ ⋅ − ∇ − = ( , ,ω, ; ) − ω ω . (3) 

Here, ) n(x,z, , E;t

r r ω is the Fourier component of the refractive index of the medium oscillating 

at frequency ω. Strictly speaking, a quantity in the frequency domain is a function of ω only, not 

of t. In the present case, however, we write n with an explicit time dependence to emphasize the 

fact that each Fourier component may slowly change with time (here, “slowly” means that the 

time scale of the variation is much longer than ω–1) as the properties of the medium evolve as a 

result of the interaction with the beam. The value of the refractive index before the arrival of the 

beam is denoted by n0. 

If the ) ( E r r r ∇ ∇ ⋅ term can be dropped, equation 3 takes the form of the inhomogeneous Helmholtz 

equation, which can be solved numerically with the “split step” Fourier scheme alluded to in the 

previous section. Roughly speaking, the split step procedure treats effects arising from the righthand side of equation 3 separately from those associated with the left-hand side. One solves the 

homogeneous Helmholtz equation over a short longitudinal step and then accounts separately for 

the effects over that step of nonlinear (and/or time-dependent linear) refraction and absorption. 

The homogenous Helmholtz equation, together with the simpler equation to which it reduces in 

the case of paraxial beams, is the subject of the remainder of this report. 

We postpone to a future report a full analysis, using the methods of Lax, Louisell, and McKnight 

(4), of the implications of dropping the ) ( E r r r ∇ ∇ ⋅ term in equation 3. For the present, we observe 

only that, since free charge is assumed to be absent, it follows from the Maxwell equation for the 

electric displacement that ) ( 2

n E r r ∇ ⋅ vanishes, or that 

(ln ) 2

1 ∇ ⋅ E = − E ⋅∇ n r r r r . 
4

2.2 Homogeneous Helmholtz Equation 

In a linear medium, the divergence of E

r

 vanishes. If, in addition, the properties of the medium 

are independent of time, then the right-hand side of equation 3 vanishes as well and equation 3 

reduces to 

( ) ( , ) 0 2 2 ∇ + k E x z = r r . (4) 

Here, k n / c = 0ω , the wave number in the medium. In an absorptive medium, n0 (and thus k) 

has a positive imaginary part. Equation 4 is the homogeneous vector Helmholtz equation. 

We designate the polarization direction of the field with the unit vector eˆ and write 

E(x,z) eˆE(x,z) r r r = . If eˆ is constant, equation 4 reduces to a scalar wave equation, 

( ) ( , ) 0 2 2 ∇ + k E x z = r , (5) 

the scalar Helmholtz equation. We now write the electric field as an envelope function times a 

plane wave propagating in the z-direction: ikz E(x,z) A(x,z) e r r = . In terms of the envelope 

function A(x,z) r , equation 5 becomes 

( 2 ) ( , ) 0 2 2 ∂ z + ik∂ z + ∇⊥ A x z = r , (6) 

where 2 ∇⊥ is the transverse Laplacian, given in Cartesian coordinates by 2 2

x y ∂ + ∂ . 

2.3 Paraxial Wave Equation 

At this point, we proceed with the slowly varying envelope approximation (also known as the 

paraxial approximation) by introducing suitable “non-dimensionalized” quantities. The problem 

itself sets a natural length scale x0 for the transverse coordinates, e.g., the size of an aperture or 

the spot radius of a focused laser beam. We use this quantity to define the scaled transverse 

coordinates: 0 X x / x r r

= . The transverse length scale x0 and the wave number k now determine 

z0, the longitudinal length scale or “Rayleigh range,”∗

 by 2

0 0 z = k x , and the corresponding scaled 

longitudinal coordinate is given by Z = z/z0. We denote by ε the ratio of the transverse and 

longitudinal length scales: ε = x0/z0=(kx0)

–1. (An alternate procedure for non-dimensionalizing 

the problem, employed by some workers (5), simply expresses all lengths as multiples of the 

wavelength λ, or equivalently, as multiples of the inverse wave number k

–1. This seemingly very 

natural choice is not the most convenient for our purposes.) 

In terms of these new dimensionless coordinates, equation 6 becomes 

( ) 2 ( , ) 0 2 2 2 ∂Z + i∂Z + ∇⊥ A X Z = r ε . (7) 

 ∗

 If one chooses x0 to be the spot radius at focus of a Gaussian beam, then z0 is the confocal parameter of the beam, equal to 

twice the Rayleigh range of the beam as conventionally defined. 
5

In equation 7, the validity of the slowly varying envelope approximation is completely manifest: 

as long as the transverse length scale x0 is much greater than the wavelength, ε << 1 and the 

approximation holds. Equation 7 then reduces to 

(2 ) ( , ) 0 2 i ∂ Z + ∇⊥ A X Z = r , (8) 

the paraxial wave equation. In terms of the original “dimensionful” coordinates (x,z) r , the 

paraxial wave equation reads 

(2 ) ( , ) 0 2 ik∂ z + ∇⊥ A x z = r . (9) 

3. A Précis of Classical Scalar Diffraction Theory 

3.1 History 

The intuitive picture presented in 1678 by Christian Huygens (6) forms the basis for the modern 

theory of diffraction. According to what is now known as the Huygens-Fresnel principle, every 

point on a wave front may be thought of as the source of spherical secondary wavelets that 

propagate through space with the same speed and frequency as the primary wave; the wave front 

at any later time is simply the envelope of the secondary wavelets. In his famous prize paper 

presented to the Paris Academy in 1818, Augustin Jean Fresnel combined Huygens’s envelope 

construction with Thomas Young’s ideas about interference. The Huygens-Fresnel principle was 

placed on a firm mathematical basis in the work of Gustav Kirchoff (7). Kirchoff’s theory, 

despite its inconsistent assumptions (8, 9) agrees extremely well with experiment in all practical 

limits. It remained only for Sommerfeld (10) and Rayleigh to remove the inconsistencies from 

the Kirchoff theory. The mathematical expression of the Huygens-Fresnel principle takes the 

form of an integral which (after a series of simplifying approximations that impose certain 

limitations on the relative size of the diffracting aperture, the distance from the aperture to the 

plane of observation, the size of the observation region on the plane, and the wavelength of the 

radiation) may be evaluated analytically in a variety of simple aperture geometries (e.g., 

rectangle, circle, and infinite one-dimensional slit). Conventional textbook treatments (11, 12, 

13, 14) typically include several of these approximate calculations. 

3.2 Rayleigh-Sommerfeld Diffraction Integral 

Although excellent discussions of the conventional integral formulation of diffraction theory are 

given elsewhere (we find presentation of (14) particularly lucid), in this section we briefly 

review the standard treatment in order to emphasize several points that surface in our subsequent 

discussion. We consider the situation at an observation point P with coordinates ) (x,z r , as 

shown in figure 1. (Throughout this discussion, we use the two-component vector x

r

 to indicate 

positions in the transverse plane.) The vector r

r

 with components (x − x′,z) r r points from an 
6

P

z

x

y

r

r

x′

r x

r

Figure 1. Diffraction by an aperture: geometry of the problem. 

infinitesimal source at coordinates ) (x′,0 r to the observation point P; the location of the source 

point ranges over the aperture. The resulting Rayleigh-Sommerfeld diffraction formula is then 

( ) ∫∫ 



 



 = ′ ′ −

plane

ikr

ik r z

r r

e E x z d x E x cos , ˆ 1 ( ,0) 21 ( , ) r 2 r r

π

, (10) 

where 2 2

r = r = z + x − x′ r r r and represents the distance between the infinitesimal source and 

the observation point, while k is the wave number, given in terms of the wavelength λ in the 

medium by k = 2π/λ. The integration in equation 10 is formally over the entire xy plane, but the 

Kirchoff boundary condition, which sets 0 E(x,0) ≡ r for points outside the aperture, restricts the 

integral to the two-dimensional aperture itself. For the case in which the aperture is illuminated 

by a normally incident plane wave, it has recently been shown that equation 10 may be rewritten 

in terms of a one-dimensional parametric integral around the perimeter of the aperture (15). 

The direction cosine in equation 10, which has the value z/r, is known as the obliquity factor or 

inclination factor. Fresnel himself recognized the need for a quantity of this kind (12), without 

which the secondary wavelets envisioned by Huygens would also generate a reverse wave 

propagating back toward the source illuminating the aperture; no such backward propagating 

wave is observed experimentally. In the Fresnel-Kirchoff diffraction formula, the obliquity 

factor is ( ) cos( ) , ˆ 1 2

1 r z + r . Interestingly enough, equation 10, including the proper form of the 

obliquity factor, can be derived entirely from Fourier transform theory, with no reference 

whatsoever to the Huygens-Fresnel principle (5). 

Equation 10 treats electromagnetic radiation as a scalar phenomenon, completely ignoring the 

fact that the various components of the electric and magnetic field vectors are coupled by 

Maxwell’s equations and thus cannot be treated independently. Within the limitations of a scalar 

theory, however, equation 10 is considered to be exact throughout the entire diffraction space. 
7

Fortunately, the limitations that a scalar theory places on the macroscopic length scale are not 

particularly severe; experiments involving the diffraction of microwaves from apertures confirm 

that the scalar theory is valid, provided the aperture size is much greater than the wavelength of 

the radiation (16), while a detailed theoretical analysis involving focused Gaussian laser beams 

shows that the scalar treatment suffices so long as the focal spot size is at least approximately as 

large as the wavelength (17). 

3.3 Preliminary Approximations 

At distances z from the aperture on the order of the wavelength, the exact Rayleigh-Sommerfeld 

diffraction formula, equation 10, can be evaluated by the method of stationary phase, which is 

described in (11). For z values of at least several wavelengths, one can take ( ) 1/r − ik ≈ −ik , in 

which case equation 10 becomes 

∫∫ ′ ′ − =

plane

ikr

r

z d x E x ei E x z 2

2 ( , ) ( ,0) r r

λ . (11) 

The highly oscillatory character of the integrand makes the direct numerical evaluation of 

equation 11 problematic; however, a recent report suggests that present-day desktop computers 

running commercially available applications software may now be equal to the task, at least for 

configurations for which the Fresnel number is not too large (18). 

3.4 Fresnel Approximation 

The standard method of attacking the integral in equation 11 relies on the expansion 

4

2

2

2 2

2

2 2

1 













 − ′ + ′ + ⋅ ′ = + − z

x x

o

z

x

z

x x

z

x

z

r r r r r r r

. 

One approximates the factor of r

2

 in the denominator of the integrand with only the zeroth-order 

term, while for the factor of r in the phase e

ikr, the three quadratic terms are retained as well. 

∫∫ 



 



 ′ ′ − ′ − =

plane

ikz

x x

z

ik d x E x

z

i e E x z

2 2

2 ( , ) ( ,0) exp r r r rλ (12) 

This is the Fresnel approximation, whose validity has been the subject of a number of 

investigations over the years (19, 20). 

3.5 Fraunhofer Limit 

If now the dimensions of the aperture are such that for every point ) (x′,0 r in the aperture, 

1

2

2

<< ′

z

k x

r

, 
8

then exp[ / ] 1 2

2 k x′ z ≈ i r and equation 12 reduces to 

( ) ∫∫ 



 



 ′ ′ − ⋅ ′ − =

plane

i

ikz

x x

z

ik k x z d x E x

z

i e E x z r r r r r ( , ) exp[ / ] ( ,0) exp 2 2

2 λ , 

the Fraunhofer approximation. 

4. Solution of the Paraxial Wave Equation 

One can easily verify by direct substitution that equation 12, the diffraction integral in the 

Fresnel approximation, solves the paraxial wave equation, equation 9, with ikz E(x,z) A(x,z) e r r = , 

thus demonstrating the equivalence of the Fresnel and paraxial approximations. It is nevertheless 

illuminating to proceed with the direct solution of equation 7 by standard methods, imposing the 

“paraxial limit” ε → 0 only when necessary. 

4.1 Helmholtz Equation as a Singular Perturbation of the Paraxial Wave Equation 

We express ) A(X , Z r in terms of its Fourier transform in the transverse coordinates. 

( ) ∫∫ − ⋅ =

plane

i X A X Z d A Z e r r r r

κ π κ (κ, ) ~ ( , ) 2 1 2

The Fourier-transformed field amplitude ) ( , ~

A κ Z r satisfies 

( ) ( , ) 0 ~ 2 2 2 2 ε ∂Z + i∂Z −κ A κ Z = r , (13) 

where κ κ

r = . Setting i Z A Z e α (κ, ) ψ (κ ) ~ r r = , we find that equation 13 is satisfied, provided that 

2 0 2 2 2 ε α + α +κ = . (14) 

There are two roots, 

( ) 2 2

2

1

2

2 2 1 1

κ ε

ε

ε κ

α = − + o − + − + =

and 

( ) 2 2

2

2 1

2

2 2

2 1 1ε κ ε

ε

ε κ

α = − + + o − − − = −

− , 

of which the second, α − , is singular in the limit ε → 0 . Returning to the physical field, 

equation 2, 
9

( )

( ) ∫∫

∫∫

= ± − +

Ε = + − +

− ⋅ −

±

− ⋅ −

plane

i X

plane

i X

e d e i Z i t c c

X Z t e d e i Z i t c c

ˆ 2 ( ) exp[ ] . .

( , , ) ˆ 2 ( ) exp[ ( ) ] . .

2

2 2 1 2 1

1 2 2

π κ ψ κ ω

π κ ψ κ ε α ω

ε

κ ε κ

κ

r r

r r

r

r r r

we see that the root α − corresponds to a solution that propagates in the negative z-direction. 

Nowhere in the region Z > 0 , either at infinity or elsewhere, is there a source from which such a 

backward propagating wave could originate, so we are justified in discarding α− as unphysical. 

Thus, whereas in the integral formulation of diffraction theory, one introduces the obliquity 

factor in order to suppress troublesome, backward propagating Huygens wavelets, here one must 

exorcise the demon with an appeal to the boundary condition at Z → +∞ ! 

In the language of differential equations, one would say that the term of order ε

2

 in equation 7 is 

a singular perturbation. Such a perturbation can completely alter the character of the solution. In 

the present case, “turning on” the perturbation introduces a second root in equation 14 where 

there was only one before and so in some sense, it is not unnatural that this new root should be 

singular. It corresponds to a rapidly oscillating solution of the differential equation, and this 

solution can be very troublesome indeed. 

4.2 Convolution Kernel for the Paraxial Wave Equation 

To obtain a solution to the paraxial wave equation, we discard the singular root and approximate 

the remaining, nonsingular root by the zeroth-order term in its expansion in powers of ε : 

2

2

1 α + ≈ − κ . 

( ) ∫∫ − − ⋅ =

plane

Z i X A X Z d e e i r r r r κ κ π κ ψ κ 2

2 ( , ) 2 ( ) 1 2 (15) 

The function ) ψ (κ

r is determined by the boundary condition on A at Z = 0 . 

( ) ∫∫ − − ⋅ ′ = = ′ ′

plane

i X A d X A X e r r r r r

κ ψ κ (κ,0) 2π ( ,0) ~ ( ) 1 2 (16) 

Substituting equation 16 into 15 and reversing the order of integration, we obtain 

∫∫ ∫∫ − ⋅ − ′ = ′ ′

plane plane

Z i X X

e e

d A X Z d X A X i ( )

2

2

2 2

2

[2 ] ( , ) ( ,0) r r r r r κ κπ

κ . (17) 

Recalling that 2 2 2 κ = κ x + κ y , we see that the integral over the κ plane in equation 17 factors into 

a product of two identical one-dimensional integrals of the form 

iZ

e G X X Z dq e e

Z

i X X

q Z iq X X i

2π 2π

1 ( , )

2

( )

( )

1

2

2

2

− ′

∞

−∞

− − ′ − ′ = = ∫ . (18) 
10

Equation 17 thus takes the form of a two-dimensional convolution, 

∫∫ = ′ ′ − ′

plane

A(X,Z) d X A(X , 0) G(X X , Z) 2 r r r r , (19) 

in which the kernel, 

( , ) ( , ) 2 ( , ) 1 1

2

2

G X X Z G Y Y Z

iZ

e G X X Z

X X

Z

i

− ′ = = − ′ − ′

− ′

π

r r

r r , (20) 

is the square of the one-dimensional convolution kernel in equation 18. Expressed in the 

“dimensionful” coordinates (x,z) r , the two-dimensional convolution kernel in equation 20 is 





 



 − ′ = − − ′

2

2

0

2 ( , ) exp x xz

ik

z

ix G x x z r r r r

λ . (21) 

Recalling that , ( , ) ( , ) ikz E x z A x z e r r = we see that equation 19 is identical to equation 12, the 

Fresnel approximation of the Rayleigh-Sommerfeld diffraction integral. 

4.3 A Remark on Convergence 

In order to ensure convergence of the integral in equation 18, 2

0 Z = z/ k x must have a negative 

imaginary part. Since z and x0 are both real, the requirement that Im[Z] < 0 is equivalent to the 

condition Im[k] = α/2 > 0, in which α is the absorption coefficient of the medium; mathematical 

convergence is thus assured in the case of propagation through an absorptive medium. Of 

course, the integral in equation 18 also converges for propagation through a lossless medium. In 

this case, we observe that a perfectly monochromatic field is a mathematical idealization. An 

actual, physical field must have been turned on at some time in the past, which one may account 

for by adding a small positive imaginary part to the frequency, so that e

–iωt

, and thus the physical 

field, vanishes as t → −∞ . This is equivalent to a positive absorption coefficient. 

4.4 Numerical Computation of the Fresnel Convolution Integral 

In practice, the integral on the right-hand side of equation 19 will generally have to be evaluated 

numerically. This is most efficiently accomplished through the use of Fourier transform methods 

using the well known fast Fourier transform (FFT) algorithm. Given a complex-valued function 

f (X )

r on the plane 2 R , we symbolically denote the two-dimensional Fourier transform of 

f (X )

r , which is a function of κ

r

, by 

 ∫∫ − − ⋅ = =

plane

i X f f d X f X e r r r r

κ κ (κ ) [2π ] ( ) ~ [ ]( ) 1 2 F

and the inverse transform of ) (

~ κ

r f , a function of X

r

, by 
11

∫∫ − − ⋅ = =

plane

i X f X f X d f e r r r r r

κ π κ (κ ) ~ ]( ) ( ) [2 ] ~

[

1 1 2 F . 

Using these relations and the well-known integral representation of the Dirac delta function in 

terms of complex exponentials, one can write equation 19 symbolically as 

( , ) 2 [ [ ]( , 0) [ ]( , )] 1 A X Z π A κ G κ Z r r r F F F − = . (22) 

Equation 22 is a prescription for calculating ) A(X,Z r : compute the Fourier transform of 

A(X, 0) r , the field at the boundary, Z = 0; multiply this by the Fourier transform of the 

convolution kernel, equation 20; compute the inverse Fourier transform of the resulting product. 

Now, the Fourier transform of the convolution kernel in equation 20 is given by 

Z i

G Z e

1 / 2

2 [ ]( , ) [2 ] κ κ π r r − − F = , 

so equation 22 is 



  = − − Z i

A X Z A e

1 / 2

2 ( , ) [ ]( , 0) κ κ r r r F F . 

This is precisely equation 15, in which ) ψ(κ

r , given by equation 16, is simply ) [ ](κ, 0 r F A , the 

Fourier transform of the boundary condition at Z = 0. Thus, we see that same procedure that we 

employed to derive the solution in equation 19 and the corresponding Fresnel kernel in equation 

20 also provides the most efficient means of evaluating the integral expression in which the 

solution is given. 

5. Solution of the Helmholtz Equation 

In this section, we derive the Rayleigh-Sommerfeld result by employing the integral transform 

methods to solve equation 5, the scalar Helmholtz equation. The calculation proceeds exactly as 

in the previous section. We take 

( ) ∫∫ − ⋅ =

plane

i x i z E x z d e e κ α π κ ψ κ r r r r

( , ) 2 ( ) 1 2 , (23) 

where ) ψ(κ

r is determined by the boundary condition on E at z = 0: 

( ) ∫∫ − − ⋅ ′ = ′ ′

plane

i x d x E x e r r r r κ ψ(κ ) 2π ( ,0) 1 2 . (24) 

The Helmholtz equation is satisfied, provided that 

2 2 α = ± k −κ . (25) 
12

The negative root corresponds to a backward propagating solution and, as before, we discard it 

as unphysical. Substituting equations 24 and 25 into 23 and reversing the order of integration, 

we obtain 

( ) ∫∫ ∫∫ − ⋅ − ′ − = ′ ′

plane plane

i x x iz k E x z d x E x d e e 2 2 2 2 2 ( ) ( , ) 2 ( ,0) κ κ π κ r r r r r

, (26) 

or, expressed in the dimensionless variables introduced in section 2.3, 

( ) ∫∫ ∫∫ − ⋅ − ′ − = ′ ′

plane plane

i X X iZ E X Z d X E X d e e 2 2 2 2 2 2 ( ) 1 / ( , ) 2 ( ,0) κ ε κ ε π κ r r r r r

. 

Utilizing the result of Lalor (21), 





















+ ∂

∂ = −

+

⋅ −

∫∫ 2 2

2

2 2

2 2

2

z x

e

z

d e e

ik z x

plane

i x iz k

r

r

r r

κ π κ κ , 

we perform the integration over κ on the right-hand side of equation 26 and so obtain 

( ) ∫∫ 



 



 = ′ ′ − −

plane

ikr

r

z ik

r r

e E x z d x E x1 ( , ) 2 ( ,0) r 1 2 r π , 

which is the Rayleigh-Sommerfeld formula, equation 10, with obliquity factor z/r. 

6. Axial Symmetry 

We now turn our attention to problems exhibiting cylindrical symmetry about the z-axis. As 

before, x0 denotes the transverse length scale characterizing the problem. We employ this length 

scale to define the dimensionless radial coordinate R = r/x0 and we use the resulting polar 

coordinates (R,θ) to designate positions in the transverse plane. 

6.1 Cylindrically Symmetric Restriction of the General Fresnel Kernel 

We impose the restriction of cylindrical symmetry on equation 19, the general solution of the 

paraxial wave equation, and on the associated convolution kernel in equation 20. Cylindrical 

symmetry requires that the field amplitude A(R,Z) be independent of the angular coordinate, so 

without loss of generality, we set θ = 0. We use the Law of Cosines to write 

2 2 2

X − X′ = R − 2RR′cosθ′ + R′ r r and so express the convolution kernel in equation 20 as 

( )



 



 ′ ′ = − ′ ′ + ′ =

2 2

0 2 cos

2 exp 2

1 ( , , , , ) R RR R

Z

i

iZ

G R R Z θ

π

θ θ θ . 
13

Inserting the above expression for the convolution kernel into equation 19 with 

d X ′ = R′dR′dθ ′ 2 , we employ the well-known integral representation of the Bessel function, 

∫ ′ = − ′

π

θ θ

π

2

0

0

cos ( ) 2

1 d e J u iu , 

to perform the integration over θ ′ and obtain the following expression for the field amplitude. 





 



 ′ ′ ′ ′ 



 



 = − ′ ∞

∫ Z

RR R dR R A R e J

Z

i

Z

i A R Z

R

Z

i

0

2

0

2

2

( ,0) 2 ( , ) exp (27) 

We write equation 27 in the following form: 

( , ) ( ,0) ( , , )

0

A R Z = dR′R′ A R′ G R R′ Z ∫

∞

, (28) 

in which the factor of R′ appears explicitly in the integrand. The associated kernel is then 





 



 ′ 



 



 + ′ ′ = Z

RR J

Z

i R R

iZ

G R R Z 0

2 2

2

( ) exp1 ( , , ) , (29) 

the analog of equation 20 for situations involving cylindrical symmetry. (By including the factor 

of R′ in the integration measure instead of the kernel, we obtain an expression for the kernel that 

is symmetric in the variables R and R′ .) 

6.2 Solution of the Paraxial Wave Equation Via Hankel Transform 

It is also possible to obtain the cylindrically symmetry kernel, equation 29, in a somewhat more 

methodical fashion, namely, by direct solution of the paraxial wave equation in two transverse 

dimensions under the assumption of cylindrical symmetry. In this case, the transverse Laplacian 

operator is given by R R R ∇ = ∂ + ∂ ⊥

2 2 1 and the paraxial wave equation 8 takes the form 

(2 ) ( , ) 0 2 1 i∂ + ∂ + ∂ R A R Z = Z R R . 

The cylindrical symmetry of the problem suggests that we express A(R,Z) in terms of its Hankel 

transform, 

∫

∞

=

0

0 ( , ) ( ) ~ A(R,Z) dκ κ A κ Z J κR . (30) 

The field amplitude in the transform domain, ( , ) ~

A κ Z , satisfies the ordinary differential equation 

( , ) 0 ~ ( )

2

2

1 i ∂ Z − κ A κ Z = . The solution is: 

Z

i

A Z A e

2

2 ( , 0) ~ ( , ) ~ κ

κ κ −= , (31) 
14

where ) ( ,0 ~

A κ is given by the inverse Hankel transform of the boundary condition at Z = 0. 

∫

∞

= ′ ′ ′ ′

0

0 ( ,0) ( ,0) ( ) ~

A κ dR R A R J κR (32) 

If we substitute equation 32 into 31, then substitute the result in equation 30 and reverse the 

order of integration, we obtain 

∫ ∫

∞ − ∞

= ′ ′ ′ ′

0

2

0 0

0

2

( , ) ( ,0) ( ) ( )

Z i

A R Z dR R A R d J R J R e κκ κ κ κ . 

Comparing this last expression with equation 28, we recognize that the κ integral is none other 

than ) G(R,R′,Z , the kernel that we seek: 





 



 ′ ′ = ′ = − ′ + ∞ −

∫ Z

RR

e J

Z

i G R R Z d J R J R e

R R

Z

i Z i

0

( ) 2

0

2

0 0

2 2 2

( , , ) ( ) ( ) κκ κ κ κ . 

The integral is given in Watson (22) and reproduces our previous result, equation 29. 

6.3 Cylindrically Symmetric Fraunhofer Kernel 

In the case of circular symmetry, the Fraunhofer condition expressed in dimensionless variables 

reads 1 /(2 ) 2 R′ Z << . In the Fraunhofer limit, the factor exp[ / ] 2

2 R Z i ′ in the integrand of 

equation 27 may be replaced by unity, which gives 





 



 ′ ′ ′ ′ 



 



 ≈ − ∫

∞

Z

RR R dR R A R J

Z

i

Z

i A R Z 0

0

2 ( ,0) 2 ( , ) exp . (33) 

Thus, the Fraunhofer kernel associated with equation 28 is given by 





 



 ′ 



 



 ′ ≈ Z

RR J

Z

iR

iZ

G R R Z 0

2

2 exp

1 ( , , ) . 

7. Example: Clipped Gaussian Beam 

We conclude with an example that is of considerable practical importance: the diffraction of a 

principal-mode Gaussian beam by a coaxial circular aperture that, for simplicity, will be assumed 

to lie in the focal plane of the beam. This problem is characterized by two radial length scales: 

x0, the radius of the aperture, and w0, the half-width of the beam at 1/e

2

 of its on-axis intensity. 

Either quantity could be used to define a dimensionless radial coordinate; we choose the former, 

setting R = r/x0. We also define W0 = w0/x0, the dimensionless quantity characterizing the extent 
15

of the beam. The clipped Gaussian beam is then described by the following Kirchoff boundary 

condition at Z = 0. 









≥

< =

−

0, 1

, 1 ( ,0)

2

0 2 /

0

R

A e R A R

R W

 (34) 

7.1 Fresnel Diffraction of a Clipped Gaussian Beam 

With the clipped Gaussian boundary condition in equation 34, equation 27 becomes 





 



 ′ = − ′ ′ − ′

∫ Z

RR

e dR R e J

Z

iA A R Z iR Z R

0

1

0

0 2 2 2

( , ) η , (35) 

where ) /(2 2 = W0 − i Z − η . After repeated partial integration and application of the well-known 

recurrence relation 

[ ( )] ( ) 1 x J x x J x

dx

d n

n

n

n = − , 

equation 35 becomes 





 



 



 



 



 



 = − − ∑

∞

= Z

R J

R

Z

Z

iR

Z

iA A R Z n

n

n

1

2

0 2

2 exp 2 ( , ) η η

η

. (36) 

Associated with the field amplitude of equation 36 is the intensity profile 

2

1

2 2 1

(0, )

( , ) 



 



 



 



 = − ∑

∞

=

−

Z

R J

R

Z

e

I Z

I R Z n

n

n

η η . (37) 

7.2 Fraunhofer Diffraction of a Clipped Gaussian Beam 

Substituting equation 34, the boundary condition for a clipped Gaussian beam, into the 

Fraunhofer result, equation 33, one still obtains equation 35, except that 2

0

− η =W , instead of 

/(2 ) 2 W0 − i Z − . We thus obtain the field amplitude of a clipped Gaussian beam in the Fraunhofer 

limit by simply setting 2

0

− η =W in our previous result, equation 36. 

7.3 Limiting Case: Top Hat Beam (Uniformly Illuminated Aperture) 

A beam that overfills the aperture ( 1 W0 >> ) provides illumination that is effectively uniform. In 

this case, ) η ≈ −i /(2Z and equation 36, the solution in the Fresnel regime, reduces to 

A R Z A e ( ) J n ( ) R Z n

n

R

i

R

Z

i

( , ) /

1

( 1) 2

0

2

∑

∞

=

− +

= . 
16

As we saw in section 7.2, the Fraunhofer limit is obtained by setting 2

0

− η =W in equation 36. In 

the present case, we have 1 W0 >> , so to obtain the Fraunhofer limit, we need only set 0 η = in 

equation 36. Doing so yields 

( )

R

J R Z A R Z iA e

R

Z

i / ( , ) 2 1

0

2

≈ − . 

The associated intensity profile in the far field is then given by 

( ) 2

1

/

2 /

(0, )

( , ) 



 



 = R Z

J R Z

I Z

I R Z , (38) 

the η → 0 limit of equation 37. Noting that R / Z = kx0r / z = kx0 tanα ≈ kx0 sinα , where α is 

the angle that the observation point makes with the z-axis as measured from the center of the 

aperture, we recognize that equation 38 is in fact the celebrated result of Airy (23). 

8. Summary 

From the Maxwell equations in a nonmagnetic, isotropic medium, one may derive a general 

nonlinear vector wave equation. If 0 ∇(∇ ⋅ E) ≡ r r r , this reduces to the inhomogeneous vector 

Helmholtz equation, which is generally written with the nonlinearity contained in the source term 

on the right-hand side. If the direction of the electric field polarization is constant, as is usually 

the case, one solves the corresponding scalar equation. In a geometry in which the transverse 

length scale is much larger than the wavelength of the radiation, the Helmholtz equation is well 

approximated by the paraxial wave equation, which is only first order in the longitudinal variable 

and thus much easier to solve. The Helmholtz equation may be regarded as a singular 

perturbation of the paraxial wave equation, and some of the difficulties arising in the solution of 

the former partial differential equation are related to this fact. 

The Huygens-Fresnel principle provides a very intuitive picture of the diffraction of light from 

an aperture. The mathematical expression of the Huygens-Fresnel principle takes the form of an 

integral over the aperture: the Rayleigh-Sommerfeld diffraction integral. The oscillatory nature 

of the integrand makes exact evaluation of the Rayleigh-Sommerfeld integral problematic, so 

one frequently employs the simplifying Fresnel approximation in evaluating the integral. We 

have demonstrated how these same results may be obtained without appeal to the HuygensFresnel principle by using standard methods to solve the appropriate partial differential 

equations: the exact Rayleigh-Sommerfeld integral may be obtained from the scalar Helmholtz 

equation in a linear medium, while the Rayleigh-Sommerfeld integral in the Fresnel 

approximation results from solution of the paraxial wave equation in a linear medium. 
17

One can employ the solutions of these linear diffraction problems to attack the more general 

problem of beam propagation in a nonlinear medium through the use of the numerical procedure 

introduced by Feit and Fleck in the late 1970s. This procedure, known as the “split step” Fourier 

method or simply as the Beam Propagation Method, separately computes the effects of 

diffraction, nonlinear refraction, and nonlinear absorption for each propagation step. 
18

9. References 

1. Perry, J. W. “Organic and Metal-Containing Reverse Saturable Absorbers for Optical 

Limiters,” in Nonlinear Optics of Organic Molecules and Polymers, H. Nalwa and S. 

Miyata, Eds. (CRC Press, Orlando, 1997), 813–840. 

2. Fleck, J. A. Jr.; Morris, J. R.; Feit, M. D. “Propagation of high-energy laser beams in the 

atmosphere,” Appl. Phys. 1976, 10, 129–160. M. D. Feit and J. A. Fleck, Jr., “Light 

propagation in graded index fibers,” Appl. Opt. 1978, 17, 3990–3998. M. D. Feit and J. A. 

Fleck, Jr., “Calculation of dispersion in graded-index multimode fibers by a propagating 

beam method,” Appl. Opt. 1979, 18, 2843–2851. M. D. Feit and J. A. Fleck, Jr., 

“Computation of mode properties in optical fiber waveguides by a propagating beam 

method,” Appl. Opt. 1980, 19, 1154–1164. M. D. Feit and J. A. Fleck, Jr., “Computation of 

mode eigenfunctions in graded index optical fibers by the propagating beam method,” Appl. 

Opt. 1980, 19, 2240–2246. M. D. Feit and J. A. Fleck, Jr., “Beam nonparaxiality, filament 

formation, and beam breakup in the self-focusing of optical beams,” J. Opt. Soc. Am B. 1988,

5(3), 633–640 (1988). 

3. Boyd, R. W. Nonlinear Optics (Academic Press, San Diego, CA, 1992), 59. 

4. Lax, M.; Louisell, W. H.; McKnight, W. B. “From Maxwell to paraxial wave optics,” Phys. 

Rev. 1975, A 11(4), 1365–1370. 

5. Harvey, J. E. “Fourier treatment of near-field scalar diffraction theory,” Am. J. Phys. 1979, 

47(11), 974–980. 

6. Huygens, C. Traité de la lumière (completed in 1678, published in Leyden, 1690). 

7. Kirchoff, G. “Zur Theorie der Lichtstrahlen,” Wiedemann Ann. 1896, 18(2), 663. 

8. Sommerfeld, A. “Mathematische Theorie der Diffraction,” Math. Ann. 1896, 47, 317. 

9. Wolf, E.; Marchand, E. W. “Comparison of the Kirchhoff and the Rayleigh-Sommerfeld 

theories of diffraction at an aperture,” J. Opt. Soc. Am. 1964, 54(5), 587–594. 

10. Sommerfeld, A. “Die Greensche Funktion der Schwingungsgleichung,” Jahresber. Deut. 

Math. 1912, 21, 309. 

11. Born, M.; Wolf, E. Principles of Optics, 6th edition (Pergamon Press, Oxford, 1980), Chap. 

8. 

12. Hecht, E. Optics, 2nd edition (Addison-Wesley, Reading, MA, 1987), Chap. 10. 

13. Siegman, A. E. Lasers (University Science Books, Sausalito, CA, 1986), Chaps. 14 and 18. 
19

14. Goodman, J. W. Introduction to Fourier Optics (McGraw-Hill, New York, 1968, reissued 

1988). 

15. Dubra, A.; Ferrari, J. A. “Diffracted field by an arbitrary aperture,” Am. J. Phys. 1999, 67(1), 

87–92. 

16. Silver, S. “Microwave aperture antennas and diffraction theory,” J. Opt. Soc. Am. 1962, 52, 

131. 

17. Chi, S.; Guo, Q. “Vector theory of self-focusing of an optical beam in Kerr media,” Opt. 

Lett. 1995, 20(15), 1598–1600. 

18. Guha, S. “Focusing by small-f-number lenses in the presence of aberations,” Opt. Lett. 2000, 

25(19), 1409–1411. 

19. Feiock, F. D. “Wave propagation in optical systems with large apertures,” J. Opt. Soc. Am.

1978, 68(4), 485–489. 

20. Southwell, W. H. “Validity of the Fresnel approximation in the near field,” Am. J. Phys.

1981, 71(1), 7–14. 

21. Lalor, E. “Conditions for validity of angular spectrum of plane waves,” J. Opt. Soc. Am.

1968, 58, 1235. 

22. Watson, G. N. A Treatise on the Theory of Bessel Functions, 2nd edition. (Cambridge 

University Press, New York, 1948), 395. 

23. Airy, G. B. Trans. Camb. Phil. Soc. 1835, 5, 283. 
20

INTENTIONALLY LEFT BLANK. 
21

Distribution

Admnstr 

Defns Techl Info Ctr 

ATTN DTIC-OCP (Electronic Copy) 

8725 John J Kingman Rd Ste 0944 

FT Belvoir VA 22060-6218 

US Army Natick Rsrch & Dev Ctr 

ATTN SSCNC-II B De Cristofano 

Natick MA 01760-5019 

US Army Tank-Automtv & Armaments 

Cmnd 

ATTN AMSTA-TR-R-MS-263 R Goedert 

Warren MI 48379-5000 

AFRL/MLPJ 

ATTN S Guha 

3005 P Stret Bldg 651 Ste 1 

Wright-Patterson aFB OH 45433-7702 

US Army Rsrch Lab 

ATTN AMSRD-ARL-CI-IS Mail & Records 

Mgmt 

ATTN AMSRD-ARL-CI-OK-T Techl Pub 

(2 copeis) 

US Army Rsrch Lab (cont’d) 

ATTN AMSRD-ARL-CI-OK-TL Techl Lib 

(2 copies) 

ATTN AMSRD-ARL-D D R Smith 

ATTN AMSRD-ARL-SE J Pellegrino 

ATTN AMSRD-ARL-SE-E H Pollehn 

(3 copies) 

ATTN AMSRD-ARL-SE-EO A Mott 

ATTN AMSRD-ARL-SE-EO G Wood 

ATTN AMSRD-ARL-SE-EO M Ferry 

ATTN AMSRD-ARL-SE-EO T Pritchett 

(8 copies) 

ATTN AMSRL-D J M Miller 

Adelphi MD 20783-1197 
22

INTENTIONALLY LEFT BLANK.