- Email: [email protected]
Calculation of forward-scattering amplitudes from a focused beam: exact values for the gaussian parameters
Calculation of Forward-Scattering Amplitudes from a Focused Beam: Exact Values for the Gaussian Parameters John Langerholc
Strassbergerstr. 10 8000 Miinchen 40 Munich, Germany
ABSTRACT The equation for radiative transfer is solved in the case of forward scattering of a laser beam that may be focused into a medium to attain its minimum in uacuo width beneath the surface. The moments of up to second order in the angular distribution of the radiance are calculated for an arbitrary point within the scattering medium and used to obtain a Gaussian fit to the radiance distribution at that point in order to reduce the mass of required computation to a small set of fundamental parameters essential for considerations of design and safety. The Gaussian parameters -u-radiance, tilt, and angular breadths of the radiance distribution-are Fourier transformed back to configuration space, where they appear as one-dimensional Hankel transforms. For calculations of the off-beam scattering amplitudes in the focal plane of the unscattered beam, a “wasp-waist” approximation is introduced, in which the radial profile of the beam shrinks to a S-function in the focal plane. A method for the improved convergence of the Hankel transforms used to evaluate the Gaussian parameters in this approximation is derived. It is then shown that these single-ray solutions, which are mainly independent of the details of the (narrow-angle) beam, can be used as Green’s functions, enabling the exact scattering amplitudes to be recovered by convolving them with the normalized profile of the unscattered beam in any plane of interest. Thus even for broad beam profiles and points near the beam axis, the computation of the exact Gaussian parameters requires only a single further integration of the relevant Green’s functions over everywhere positive, i.e., nonoscillatory, functions. The asymptotic behavior of the Green’s function for the irradiance is studied at small radial distances r, and an approximative formula for this function is derived.
1.
INTRODUCTION
It is well known that atmospheric scattering degrades the performance of laser systems, reducing visibility, attenuating the direct rays, and blurAZ’Z’LiED MATHEMATICS AM) COMPUTATZON 48:101-137 0 Elsevier Science Publishing Co., Inc., 1992 655 Avenue of the Americas, New York, NY 10010
101
(1992)
0096-3003/92/$05.00
102
JOHN LANGERHOLC
ring contours. In atmospheric applications of high-energy laser beams, scattering can not only reduce the amount of power delivered to its goal, but can even pose a safety hazard by rerouting significant amounts of energy outside the intended path of the beam. This is true not only for observers situated near the axis in the path of the broadened beam and looking in the direction of the transmitter, but also for observers at the transmitter site looking down the length of the beam. The relevant effects of back scattering from a focused beam will be studied in a follow-up work; the foundations are laid here in the calculation of forward scattering amplitudes at arbitrary depths. Although the solution of the transport equation in the small-angle approximation has been known for some time [1-lO],previous investigators have experienced numerous difficulties in numerical evaluation of the inverse Fourier/Hankel transforms involved in off-axis calculations of forward scattering amplitudes. Various methods have been suggested for evaluating the fourfold inverse Fourier transform required to determine the angular distribution Z(&r, z) for a given point in space (r, z). A Taylor series expansion of the radiance Z(&‘,r, z) about c$’ = 4 in the convolution term has been proposed by Fante [3]. Significantly different results were obtained by Tam and Zardecki [6] by expanding the exponential en in a power series and integrating term by term. Even in the relatively simple case of a conical beam diverging with Gaussian angular spread from zero width at the origin, however, the computational effort required for the evaluation of the nfold integral occurring in the nth term proves to increase so rapidly that the latter authors break off the summation when the last summand becomes less than 5% of the total attained. This results in taking from 3 to 6 terms, depending on the axial distance. The resultant expressions for the nth term are sufficiently complicated even in the calculation of the irradiance, which no longer depends on the angle 4, that special methods have been devised to improve convergence [ll, 121. In [II] ,a resummation of the series was achieved, leading to an integral over a K, function with exponential falloff, but complex argument. In [12],summation was carried out by a rather complicated computer program. In the light of the present, much simpler results, neither of these procedures appears very satisfactory. 1.1. OjfBeam Zrradiance In addition to the difficulties already experienced for a conically diverging beam, the desire to calculate in the waist plane of a convergent, sharply focused exacerbates the problem touched on in [9] ,where
in such calculations off-beam amplitudes laser beam further the kernel r~Ja(r~r)
Forward Scattering Amplitudes
103
oscillates with increasing amplitude and is only very slowly damped in the integration variable 7 for radial distances r larger than the beam width. 1.2.
Data Compression
To calculate at a given point the radiance, which depends on both the position of the point and orientation, it is theoretically necessary to evaluate a four-dimensional inverse Fourier transform. For many considerations, including those of laser safety, where it is a question of maximal exposure, it is not necessary to calculate the exact functional dependence of the radiance on the orientation angle in those cases where the breadth of this distribution is so small as to render all the light scattered to a given point effective. This condition is expected to be met in the case of narrow-angle forward scattering, and in the case of sharply peaked back scattering to be treated in a subsequent work. For this reason, a theoretical approach will be followed here which first reduces the computational volume from a fourfold inverse Fourier transform at each point (r, z) and every orientation 4, to a set of twodimensional transforms, condensing the essential information into a set of relatively few parameters basic to the angular distribution, including irradiance, average orientation and angular breadths in the radial and tangential directions. Since these measures figure as parameters in a Gaussian representation of a beam, they are referred to here as the Gaussian parameters of the scattered radiance distribution. In cases where these parameters in themselves are not sufficient, the radiance can still be estimated as a function of the orientation angle Q, by calculation from the Gaussian distribution determined by the calculated parameters, as discussed in the following subsection. 1.3. Focused Beam The configuration to be studied is illustrated in Figure 1, showing a cylindrically symmetric laser beam focused into a scattering medium, where it attains a minimum width of Ri, = l/y at a penetration depth z = zr and then diverges again at greater depths. Although originally intended for this special case, the calculations are general in the sense that zr can be chosen negative or zero to place the focus outside the medium or at its boundary. The latter case is the one that has been studied almost exclusively up to now. The scattering phase function will be assumed as in [5-g] to have a Gaussian forward peak; residual isotropic scattering will be treated in subsequent works by means of a two-component scattering theory [lo].
JOHN LANGERHOLC
104
-----
_-._,-.- .._.-.- .--. -.-.
Input
FOCd PlalW
Plane LO
z=z, FIG. 1.
2.
RADIATIVE
TRANSFER
THEORY
The basis of the present calculations will be the radiative transfer or transport equation derived according to geometrical optics from energetic considerations. Interference and polarization effects are neglected; the fundamental quantity to be solved for is the radiance or specific intensity distribution Z(+,x). The solution of the basic equation in the small angle approximation will be reviewed to assist in generalizing the known solutions to the case of a beam not necessarily focused at the input plane. The Transfer Equation The basic equation for the radiance derived from consideration of a differential cylindrical volume element with generator in the direction of the unit vector 1, is [l-4] 2.1.
l,*v,q~,x)+ uJ(4,x) = Uf// +wy(4’,x)
dfl,,.
(2.1)
Here P(+, 4’) is the ph ase function, giving the probability that a ray will be scattered from the angular direction 4’ to q5. It is normalized to unity: s 4a
P(&,&)
dn,
= 1;
(2.2)
and the scattering coefficient is designated by ur, the amount of power diverted out of a differential volume dV = dA* dz being given by = ufzdv= dZSCatt
(ufdz).(Z-dA).
(2.3)
105
Forward Scattering Amplitudes
Similarly, the extinction coefficient gives the total power lost from the beam when traversing the elementary volume: dP’=
a,ldV.
(2.4
In the multicomponent scattering theory of [lO],it is composed of the various scattering coefficients corresponding to forward peak, isotropic component, and rear peak, along with the absorption coefficient: u,=uf+a,+u,+u,.
(2.5)
For the present work dealing with forward scattering, two terms in this equation will be relevant.
only the leftmost
2.2. The Small Angle Approximation General solutions for the radiative transfer equation are difficult to find; and it is customary in the case of forward scattering of a collimated laser beam to make the following approximations, based on the assumption that the bulk of the scattered light propagates in directions deviating by only a small angle from the principal axis 1 z. of the laser beam [l-9]. The coordinates are represented as x = (r, a);
and the angle C#is spanned by a two-dimensional
(2.6) vector
1,=(&l),
P-7)
where
The angular integration thus turns into a planar integral:
In the remainder of this paper, the limits of integrals that are taken from - 00 to 00 will be omitted. Integrals without limits should be understood as such.
106
JOHN LANGERHOLC
The phase function is also considered to depend only on the magnitude of the vectorial difference of incident and scattering angles:
p(w)= With these approximations becomes
=U
f
q M-4’1).
(2.9)
and changes of notation, the transfer equation
JJ
P(+ - ti’)I(+‘,r,
(2.10)
2) d24’.
Fowi.er Transformation of the Transfer Equation 2.3. In this form, the natural method of solving the equation is to perform a fourfold Fourier transform in the coordinates 4 and r leaving the axial coordinate z as an independent variable:
~^(t,q, z) = J
J
J
J
(2.11)
d2~d211e’(0’~+II”)Z(~,r,z).
The difficulty of solving the equation in physical space is thereby transformed into the difficulty of evaluating the fourfold inverse Fourier transform to obtain the radiance distribution:
z(4,rF 2)=(27r-4 J J J J
d2~d2rle-‘(E-~+?‘r)i(f,tl,
z).
(2.12)
It turns out not to be advantageous to transform the equation in all variables at once, since the factor C$ would turn into the differential operator - i V{. 2.4. Fourier Transformation of Geometrical @tics The basic idea of the solution is to transform the ray propagation condition of free space Ar = +Az into a substitution for the conjugate variables which will allow the solution to be extended from the input plane z = 0 by a simple variable substitution. The equivalent of the ray propagation condition will turn out to be the substitution
E+f+zv
(2.13)
107
Forward Scattering Amp&&s
inserted into the transformed radiance at z = 0. This will not only allow the solution to be written down relatively simply from the transformed boundary values but will also provide insight into the problem of representing, and solving the transport equation for, a focused beam. 2.5. Spatial Fourier Transformation To put this statement on a rigorous basis, the transport equation is transformed at first in the two radial coordinates only, giving [4]
(2.14) The constant factors of the unknown on the left-hand side can be converted into an exponential factor by introducing the intermediate unknown I’, which is defined by:
i(4,v,r)=exp(- o,z+i+-rlz)i’(+,q,z)
(2.15)
This simplifies the equation to
13,i’(c$,q, 2) =
uf//
P(+- q5’)e-i(“-“‘)‘qZ
x i’($‘,q, z) &jJJ. 2.6.
Solution of the Transport
(2.16)
Equation
Performing now the Fourier transformation in the angular variable and using the convolution theorem, we find
where (2.18) The solution of the last differential equation can be immediately written down; it is:
Q,% z) = e”(f*-~*t)?(f,tj,O),
(2.19)
108
JOHN LANGERHOLC
where
and
qt,rl>q = @,%O). From
the above definition
of I’, it follows that
=exp(-o,z)i’([+qz,q,z). Inserting
the solution obtained
(2.20)
for Z’ yields
i(l,g,z)=exp(-~~z+62(S+11z,-rl,z))i’(~+?z,r1,0). Using the equality
(2.21)
of Z and I’ for z = 0 and the fact that
n(o,?z,-~r,~)=u~J’~(E+a(z-z’))dz’ 0
=(J
f
J
‘p(E+Ilz’)dz’=n(E,~,z),
(2.22)
0
we have, finally, i($,~.z)=exp(-
a,~ + Q(E>l),
z))Z(f
+ MLO).
(2.23)
This may be expressed in words as follows. On propagating through a depth increment z, the Fourier transform of a given input radiance distribution is changed in three ways. The conjugate variable E of C$ is replaced by t: + gz; and the entire expression is multiplied by the extinction factor and by e”.
Forward Scattering Amplitudes
109
Contained within this general prescription for translating the Fourier transformed radiance in the z direction is a substantiation of the claim that the radiance of a beam propagating in uacw, can be calculated from its Fourier transform at z = 0 by replacing the conjugate angular variable [ by E + Qz. It will be used in the following discussion to derive a representation for a beam focused into a medium. This behavior is also reflected in the argument E + 7 z’ occurring in the Fourier transform of the phase function. 3. 3.1.
THE
PHASE
Representation
FUNCTION of the Forward
Peak
The conventional representation of a forward peak in the phase function is a normalized Gaussian in the transverse components of the (tangent 04 the angle [S-9]:
P( 4) = :exp(
- (~‘4’)
JJ
P(+)d2$=1.
(34
The Fourier transform is also Gaussian: F(E)
=exp( - t2/4a2).
(3.2)
This representation of P leads to the following expression for the exponent of the scattering function:
Q(e,v,~)=u~~‘exp(which is expressible
(f+qz’(‘/4cY2)dz’,
(3.3)
in terms of standard functions:
X (erf(
nz+~~se)
-ej
F)),
(3.4)
110
JOHN LANGERHOLC
where 0 is the angle between the vectors t and q:
P-5)
t’q=[qcose. For [ = 0, this simplifies to
(3.6) 3.2. Multi-Gauss Peak In the case of certain phase functions, it is desirable to take a multiple forward peak into account [IS-141 ,in which a multi-Gaussian phase function is constructed as a convex linear combination of single Gaussian peaks of different widths: P(f#l) = r-l
5
exp (-
K,c$
m=l
c&“)>
(3.7)
where
5
K,,,=l.
m=l
The result on the solution for the scattering amplitude is that Q is replaced by the same linear combination of the Q functions calculated for the individual scattering widths:
P-8) 4.
REPRESENTATION
OF THE
INPUT
BEAM
To find the forward scattered radiance from a laser beam directed into the medium, it is necessary to insert the Fourier transformed radiance of the beam at the input plane z = 0. Gaussian Beam The radiance distribution for the laser beam used in [9] is taken to be Gaussian in both angular and radial variables in the input plane
4.1.
Z(f$,r,O)
= FP2y2ap2exp(
-
P2+2- y2r2),
(4.1)
111
Forward Scattering Amplitudes
where F is the total input power of the beam. The Fourier transform is again Gaussian: i([,q,O)
= Fexp(
-
52/(2fl)2-71~/(2’~)‘]-
(4.2)
According to the general solution found above for the radiative transfer equation, the solution propagates in vacua according to the substitution
i.e., in the absence of scattering and absorption, the radiance in the plane at depth z is given by I^(~,ll,~)=Fexp(-(f+112)~/(2p)~-8~/(27)~).
(4.4)
4.2. Inward-Focused Beam The input beam described in the preceding section has its minimum breadth W(0) = v&/y in the input plane z = 0 [3, lo]; it spreads with increasing z, having a width at any given depth of
R(z)”= iw”( z) = l/y2
+ ( z/q2.
(4.5)
This solution is symmetrical about the input plane z = 0; it represents a beam converging to a relative focus at z = 0 and then diverging at the same angle l//3. To represent an unattenuated beam focused at a given depth zr beneath the surface, it is merely necessary to translate this unscattered beam along its own axis by substituting z - zf for z:
P-(&q,
z) =
Fexp{ - (t + II( 2 - ~f))~/(2/3)“gels)-
(4.6)
Thus the input radiance of the infocused beam is
Pm(t&o)
= Fexp{
- (E - Irzf)“l(2b)“-
q2/(2y)“].
(4.7)
Its width at the input plane is given by
R(0)”= iw”(o) = l/y2 + ( zf/q2.
(4.8)
JOHN LANGERHOLC
112
By means of these relations the distribution can be fully determined by the widths in the focal and input planes and the separation zf between the two planes: r=JUw(z,),
p=
v%Zr/J{W2(0)-w”(q)} . from a
Radiance Scattered
4.3.
(4.9)
We
are now in a position Fourier space for the scattered the input
plane
as boundary
(4.10)
Focused Beam
to write down the radiance by inserting value
in the general
complete solution the beam radiance solution
found
in at
above.
The result is
f(l,q, z) = Fexp(Xexp{ 4.4.
Q(~,T, z))
uez + -(I
+ v( z - 9))“/(2S)‘-
v~/(~Y)‘}.
(5.11)
Representation of the Beam in Configuration Space
The inverse Fourier transform of the unattenuated beam at a given depth z, measured again from the focal plane, can be performed analytically after the variable
substitution:
d2T = d2f
7=c+7jz; which
acts in the exponential zbeam(‘pr,
z)
=
to separate
the integration
(2T)-4F/- 1 dz,,/ J X exp( - r2/4P2
= Fb2y2T-2
(4.12) variables:
d27e-‘(‘.*+‘.-““.&)
- q2/4y2)
exp( - (/3+)‘-
With r and z constant, the &dependence Gaussian form by completing the square:
y2(r-
2#)‘}.
may be explicitly
(4.13) rewritten
in
~82+(~y)2)~2-2zy2r~d=(P”+(~y)2){~~y2rl(B2+(w)2]]2
- y4z2r2/(
P2 + ( z~)~].
(4.14)
Forward Scattering Amphuh Insertion
113
above leads to the expression
where
Zy)2)= l/~2’Y2R2( z,
(A+)” = l/( f12+(
zry2=
zry2(A+)‘,
(4,) =
P2 + ( zY)2
R”( 2) ‘iW2(
z) = ( z//3)2+l/y2.
This may be taken
as a motivation
dependence
radiance
of the
for the Gaussian
distribution
fits to the angular
to be developed
in the
next
section.
5.
GAUSSIAN Rather
than
FIT
TO THE
producing
ANGULAR
tables
DISTRIBUTION
of numerical
values
by means
of an
intensive computational procedure that evaluates the radiance in every angular direction + for a given point (r, z), the orientation, magnitude, and widths of the bell-shaped angular distribution will be calculated given point. This will allow the functional form to be modeled equivalent two-dimensional Gaussian distribution of the form
=
-
W)A4
?rA+
r
.9
exp( - (4,
-(+,))“/A@
- +82/A&),
at the as an
(5.1)
and enable an adequate fit to the distribution. The physical significance of the parameters is as follows: Z(r, z) is the irradiance, the angular integral of the radiance at the given point. The mean or expectation value ($,) of the radial component of the orientation angle 4,. is a measure of the tilt of the radiance lobe outward along 1, away from the principal direction of propagation of the laser beam 1 t. The
114
JOHN LANGERHOLC
corresponding peripheral component of the tilt vanishes by symmetry. The two widths of the distribution A+, and A+, are analogous to the inverses of the parameters CY,P, and y for the phase function and the input beam. The advantage of such a Gaussian representation is that the parameters are independent of the orientation angle 4; hence only two-dimensional Fourier transforms need to be evaluated, rather than the numerically more demanding four-dimensional transforms of the radiance function. 5.1. Moments of the Angular Distribution The fundamental parameters of the fit can be calculated from the moments M m,“[ I] up to second order (m + n < 2) in the angular variables, where Mm,,[f]
-/
/
(~~)“(~y)“f(~)d2~.
(5.2)
To avoid unnecessary notational complications, we take advantage of the cylindrical symmetry of the problem and align r along the r-axis by setting y = 0. In so doing, 4, becomes 4, and & becomes c$,,. The zeroth order moment is simply the irradiance:
M,,,[ I] = /
1 Z(hr,
It is obtained by inverse Fourier
Z(r, z) = (27rP2/
z) d24 = Z(r, 2).
(5.3)
transformation / e-““‘Z(q,
z) d2q
(5.4)
of the solution obtained above
i(% 2) = i(0,rl.z). The expectation values of the angles are obtained corresponding first moments by the zeroth moment:
(4%) = M,,o[11/%lJ[ 11. As remarked previously, (&J will vanish by symmetry.
(5.5) by dividing the
(5.6)
115
Forward Scattering Amplitudes
If the Gaussian fit I, is inserted for I in the calculation of the second moment, we have
x exp( - (4, -
(~J,))~/W~ - 4,‘/&,‘] I(r, 2)
= { (4,)”++br)‘}Mo,o[~] . Thus the condition for matching the second radial moment, M,.,t
(5.7) M,,,[ &] =
Zl,yields +4r)”
=
M,,o[
11
Po,o[
I]- (4,)“.
(5.8)
Similarly, but more simply,
,(A40)2= ~0,2[~1/~0,0[~1-
(5.9)
5.2. Moments in the Fourier Representation If the Fourier transform of the distribution is known analytically, the moments are readily calculated by substituting - iV[ for 4 and evaluating the appropriate derivative of the distribution at f = 0, which corresponds to integrating over all values of 4. Hence the moments of interest, still Fourier transformed in the spatial coordiante r, will be given by M,,,[ i] = - ial(fJ,tl, M,,,[
~)/a,$,,
i] = - d21(0, 1, ~)/a[,“>
M,,,[ i] = - a21(O>~, ~)/a$
(5.10) (5.11) (5.12)
To obtain these quantities as functions of the spatial coordinates, it suffices to perform the two-dimensional inverse transformation in 7. 5.3. Di&mntiating an Exponential The relevant part of i for the differentiation in E can be represented as ef(z), where f(E)S-{t+7(z-
2/)}2/482+n(t,a,z)+Const(r,z),
(5.13)
116
JOHN LANGERHOLC
and where the last term does not depend on 1. Because of the exponential form, we may write aefcQ/i3E, =
fx(t)ef([),
(5.14)
t )) efl%
Peflf)/8Ei
= { fX( f)2+
f,, J
Cr2ef(l)/a[i
= ( fY( e)‘+
f,, &f
))@).
(5.15) (5.16)
To obtain the moments as functions of the spatial point (r, z), it will then suffke to insert the appropriate differential expressions into the Fourier inverse transform integral for the irradiance, the correspondence being:
$-$I, 1, 2)
M 1.0
*
M 2.0
--f*(o,l,2)2-f,,.(9rl,z)
M 0.2
--f,(o,r,,r)“-~~,,(3,,,r)
-
(5.17)
With the expression given above for f, we find f,(e)=+,++-
~~)}/2P~+an(~,rl,z)/aE.,
(5.16)
where
At f = 0 the last integral can be evaluated analytically:
aqo,tl, 2) = - 7 at,
J,‘exp{
-(qz’/2o)“)
d(gz’/2a)2
(5.19)
117
Forward Scattering Amplitudes
For the first derivative, then, we have ,(,)JJz-zf) x
-~[l-e”p(-(~z/201)2)).
2P2
(5.20)
The second derivative is f,,,(E)
= -1/2P2-
Q(k%
z)/202
+ qzz’)2exp[
- (4 + ~z’)~/~cY~} dz’.
(5.21)
The last term can be simplified for ,$ = 0 by partial integration of
J‘(
qx
z’)‘exp(
-
( 1Z’/2~)2]
dz’
0
= -
exp( - (qz’/2o)‘)
di-
zexp[ - (qz/2a)2]
so that
- gexp{
(5.23)
- (qz/2cu)2).
Similar results are obtained for the derivatives with respect interchanging qX and vy. 6.
GAUSSIAN FIT PARAMETERS
AS HANKEL
to t,
by
TRANSFORMS
The quantities to be calculated for a given vectorial separation r from the beam axis at the penetration depth z have been derived above as
118
JOHN LANGERHOLC
functions of the Fourier following form:
conjugate
variable
q, and can be written in the
where i is the Fourier transform of the irradiance i(‘l.z)=F.expI-Sr1’/R2(z)-u~~+n(?,~)},
(6.2)
and R(z) is the beam width at depth z: R”( 2) = (2 - z_J2/4p2 + l/y?
(6.3)
The transformed n-radiance depends only on the magnitude of q, whereas the f(q, z) may contain combinations of r,rr, 7: and $. The asymmetry arises because r = (x, 0) is chosen parallel to the r-axis so that r = X. The quantities Q, including irradiance, tilt, and angular breadths of the radiance distribution, are calculated as inverse Fourier transforms of the Q:
Q(r? 2) = (2~)-“/
/
d”?e-““.‘@~,
z),
(6.4)
where
d2q = qdvde. In the sequel, for the sake of simplicity, the calculations will be performed in the focal plane z = zf. Hence the depth coordinate will be understood to be zf and suppressed in the functional arguments of Q. Cylindrically Symmetrical f If Q is independent of 8, the angular integral becomes a Bessel function of zeroth order: 6.1.
I
0
2*
de e-iqr.cose = 2rJo(53-),
(6.5)
119
Forward Scattering Amplitudes
hence
Angular Dependent f If &II)= (- irlx /~)“drl), then
6.2.
JJd211e-i1).'6(1)=(d/dx) J
= p/w
J
d2ve_iq’rg-“q(q)
J J 7l~a9(~)a-“Jo(v)
and hence
If, fin&,
rewritten
J
= 9(v) sin28, then the O-integral can be
O(v) = (71,/9)~9(9)
2* e
-'*"OSe,in2ede=(i*r)lJo26ine
&-ivcos0
0
= (i/qr)
de=2?rJ,(vr)/qr, J2=e-isrcoseCOSe
(6.6)
0
and hence Q(r)
J J d2~e-iq.rO(~)
= (2r)-"
=-
1
2*?-
=& Joo~~~++7(+-) [QJI Jmds./l(sr)s(v) 0
120
JOHN LANGERHOLC
6.3. Fallofl Properties of the Q(v) Difficulties in calculating the Jo-transform for the irradiance were encountered by previous researchers [5-91 because of the extremely slow falloff of the function &q/r) = i(v/r, z) in the variable 7. Combined with the factor v arising from the measure d2~ = r] dv de and the only very slightly damped oscillations of the Bessel function, this means that a direct evaluation of the integral
P-7) as it stands will involve integration over a large number of lobes of alternating sign, which initially increase rapidly in absolute value (as Jq). The method is recommended in [9] only for radial distances from the beam axis on the order of the radius of the unscattered beam or less, since the factor exp( - +q2R2( z)) will provide sufficient falloff in these cases to limit the required range of integration, hence the growth of qJa(s) within the range of interest. 6.4. W-Beam Zrradiance For such applications as calculations in the focal plane, this state of affairs is quite problematical, however, since the irradiance and other parameters of the Gaussian fit are to be evaluated at radial distances for which the unscattered beam has practically zero amplitude. After substitution of q/r for q, the falloff contribution from the beam width becomes which is practically ineffective for r] of moderate exp( - $R2( z)/4r2), size. The situation becomes dramatic if at the desired depth zf, the focused beam attains a minimum unscattered radius R( z,-) = l/y on the order of 1 cm, for instance, whereas a typical radial distance of interest would be on the order of r = 10 m. Falloff from the damping factor begins to become appreciable only above q = 2 r/R(O) = 2000. 6.5. Excision of the Unscattered Beam Since its value is not only known, but in this special case is expected to be practically zero, the contribution from the direct beam need not be calculated along with the scattering amplitude, but may be subtracted from the Fourier transform of the irradiance before evaluating without materially changing the value of the result to be computed. This value is obtained by setting Us = 0, hence Q = 0, in the above expression for i, while leaving the original extinction coeffkient a, unchanged.
121
Forward Scattering Amplitudes In the remainder
of this calculation,
tion to the irradiance
thus, only the scattering
contribu-
transform
il(~,z)=F.exp(-~R2(a)r12-u~z)(e”(?.Z)-1)
(6.8)
will be transformed back to configuration space. Since qQ(q, z)+ aa,& as ~+oo, the last factor on the right hand side will fall off as l/q. After the substitution
of q/r
for 7, the exponent
becomes
n(q/r, 2) = (ckwfrJ;;/7j)erf(7pz/2ar). The falloff due to the factor (e” - 1) can be estimated, following inequality: lea-11
(6.9) if desired,
from the
***)
(6.10)
The Wasp- Waist Approximation
6.6.
In view of the numerical facts just discussed, the factor exp( - aR”( zr)v2), which provided the sole damping in the original expression, can now be completely ignored, since within the effective range of the Bessel function, which falls off as l/~/q, its value will be essentially equal to 1. The convergence of the integral is insured by the fact that the maximal values of lo oscillate and decrease monotonely in value to zero. The result after subtraction
is invariant with respect
to the details of the
beam $(r]/r, as could be expected
Z) 2: F*exp(
- u,z)(e”(““XZ)-l),
from the small angle approximation.
(6.11) These
details
will not figure in the angular integrated irradiance but, by virtue of their appearance in the f(q/ r ), will enter into the results for the tilt and angular breadth of the radiance distribution. Geometrically, this approximation is equivalent to assuming a conical beam with radius R( zr) = 0 in the focal plane. For a sufficient radial distance within this plane, the constriction of the beam will have practically no effect on the amplitude of forward scattering at distances large with respect to its focal width. (The calculated value at r = 0 will nevertheless diverge, clearly showing the limit of this approximation.)
122
7.
JOHN LANGERHOLC
RECOVERY
OF THE
EXACT
IRRADIANCE
The usefulness of these calculations is not, however, limited to the focal plane or even to calculations of an inward focused beam. On closer inspection, it can be seen that it need not even remain at the level of an approximation. In the sense of function theory, the omitted factor exp( - aR”( 2)~~) represents a convolution of the inverse transformation of the single ray factor (en - 1) with the normalized unscattered beam profile in the plane at depth z. For off-beam values of r, where this function varies slowly with respect to the radius of the beam profile, it represents only a negligible averaging (or smearing in the sense of distributions) of the Green’s function from the single-ray solution and thus has a negligible effect. For calculations near the axis or at penetration depths for which the beam radius is appreciable, the error made in the approximation can be undone by subsequent convolution of the single-ray solution lr in configuration space with the reduced beam intensity normalized to unit power flow:
II(”
JII(r’.ZIR=O)IG(r-r’,z)dZr’;
(7.1)
where
&(r, z) = Z(r, z 1uf = O)/F( 2) = exp( - r2/R2( Z))/7rR2( z),
IS
(7.2)
Zb(r, 2) d2r = 1.
This procedure is illustrated in Figure 2. Here again it is understood that to obtain the reduced beam intensity I,, the extinction coefficient ue is maintained at its original value; the scattering coefficient ur is set = 0 as an independent variable. In the case of a Gaussian beam profile, the angular integration can be performed analytically. It leads to a convolution of a product of a Gaussian
123
Forward Scattering Amplitudes
FIG. 2.
with the Bessel function of imaginary argument I,:
(‘-3) where
i( r, I-‘) = f
s
* dtl Zb(r-r’, --*
z)
dBexp{2rr’cosO/R2(z)}
(7.4)
124
JOHN LANGERHOLC
Once the Green’s function Zf(r, z 1 R = 0) is calculated for all values of r, the total scattering amplitude for any axial distance from the beam axis is obtained by a single one-dimensional integral over a positive integrand. The l/r singularity of Z to be investigated in more detail in the next subsection is neutralize d by the measure dr” = 2r’dr’, ensuring an everywhere finite integrand. 8.
IMPROVEMENT OF THE GREEN’S FUNCTION
CONVERGENCE
FOR
THE
In spite of the fact that the above elimination of the direct beam now makes the integrals [Q. 1- .3] convergent even without the damping of the beam profile, the speed of convergence is nevertheless still unsatisfactory for numerical evaluation. In the integral expression [Q. 11, for example,_the factor of r~remaining from the measure d2v cancels the falloff of Q = Zr at T,J+~ and leaves only the 7 - ‘/’ falloff of the oscillatory Jo. The wasp-waist approximation, however, which eliminates the slow exponential damping, permits the following rewriting of the Hankel transforms, using Bessel’s equation and integration by parts, to increase the speed of convergence and thus make a rapid numerical evaluation feasible. 8.1. Zmproved Convergence of the Jo-Transform From the Bessel equation of order zero
d2Jo( rl) + 1 dJ,(rl) + Jo{11)= 0, h2
rl
(8-l)
ds
we may write
J&l) = -
d%( rl) 1 dJo(v) 7
-
-
11
Setting o( q/ r) = f(v) and inserting Equation 27rr2Q(r)
= - /j+)
d&(rl)-
-
drl -
(8.2) \ I
.
[8.2] into [Q.l]
somf(
yields
v)dJo( rl).
(8.3)
A partial integration converts the first term on the right into
-I;(q)qf(q)lo+~m;l(q){f(n)+?f.(q (8.4)
125
Forward Scattering Amplitudes
The first term vanishes as long as f(0) is finite and rlf(q) remains bounded as 7 + 03. The first term of the integral cancels the last integral of the equation above, leaving
27rr2Q(
4 = /j.f’(
7) 4o( 7)
= lo(4)vf+)lrn-
0
/“lo(v)
+Lf(?)).
0
P-5)
Assuming rlfl(q) is bounded as well, as will turn out to be the case in the forthcoming discussion, the first term again vanishes, leaving
t 0.1’1
8.2. Calculation of the b-radiance To calculate the scattered irradiance to insert into [Q.l’j is
I#-, zr), the appropriate function
(8.6) where C= F*exp( - u,zf),
w(q)=fi(7j/r,zf)=
VG arui
(8.7)
erf(s).
Now
f(v) = C*exp{w(s)}w’(rl)
(8.8)
126
JOHN LANGERHOLC
and
J,‘;fexp( - (~z’,~cx~)~) dz’
w’( 7) = @f
= uf zfexp ( - (w/242)
(8.9)
- o( rl)
so that
X
(-
(~/24~))
(8.10)
Hence the Hankel transform for the scattered irradiance is
X
-qw’(q)“+~‘(q)+$exp-(qn/2ar)~
dq.
(8.11)
8.3. Asymptotic Behavior of the k-radiance at Small r The most rapidly increasing contribution to the scattering n-radiance Zf(r, z) as r + 0 comes from the second term in the exponential expansion: ‘;( v,z)
= C( en(q~a)- 1) = Cj$
Q(v, z)m/m!,
with C= F*exp( - u,z).
(8.12)
Forward Scattering Amplitudes Inserting only the expansion,
127
m = 1 term
of this series into the inverse Bessel
=(oiC12~)SIdz’Smexp(-(llt’/2ar)‘)rllo(rg)d~. 0
(8.13)
0
and exchanging orders as has been done, we may evaluate the q-integral as a Laplace transform after a substitution of variables t =
?j2
p = ( 2’/2 Cq2 a=
r2/4
giving
00 =
J0 =
exp( - pt)Jo(2Jat)
(2ol/z’)‘exp(
dt =
p-le-‘lp
- (W/T?)“),
(8.14)
hence
I,(r, 2) = (Cu__cx2/7r)J,‘dz’exp(
-(rcy/z’)2)/z’2.
= (Cu__cz/7rr)[r,zdxexp(
= sexp(
- u,z)erfc(ar/z).
- x2)
(8.15)
JOHN LANGERHOLC
128
Thus the single-ray Green’s function has an asymptotic behavior proportional to 1 /T near the axis. 8.4. Calculation of the Tilt The first moment of the angular distribution
is calculated
by inserting
i)(S,Zf)=-if,(o,q,Z)i(‘1,Zf(Y=03),
(8.16)
where
f,(OA Zf) = - y-
(1-exp(
-(~.z@Y)~))>
(8.17)
into [Q.Z] with n = 1: M,,,[ I] (r, Zf) = bw(
?-, 4
=- J 1
27rr2
=&7(a/r)--;i?O
aJ( 7)
(8.18)
where
9(rl)=~(n)/(-iq,/q)=-~{1-rxp( -(~)“))f(vpzf) (8.19) 9(rllr)
Fuf r = - -exp(rl
2
a,zf)
e4rl)
.ii)
.
Inserting above yields
(Wf
with
(r, Zf.) =
(8.20)
129
Forward Scattering Amplitudes
Inserting Jl(s) dv = - 4drl) and integrating by parts improves the
convergence
of the expression:
x
d(n)(l--P{ - (Wf/2ar)2))
{
vf2 +-
2a2r2
exp
1( -
Wfl
2
ml-
)I) 2
(8.22)
8.5. Alternate Expression for ltilt If instead the limiting value 1 is subtracted from the exponential e“ to improve convergence, then there will be an extra term to Itilt, obtained by setting e* = 1, and equal to
J,-{1 -exp(
- (lz~/2ur)2]}lr(~)
ds
(8.23)
From standard tables of Hankel transforms, the last integral is found to have the value l-exp( where a=(+-/cxr)
‘tilt(
-l/a),
2. Hence the alternate expression for evaluating Itilt is
r,=exp(- (arlzf)2) +l-(I-exp(
-(‘lzf/2orr)‘))(eO(*)-I)J,(rl)dtl
(8.24)
An interesting feature of this approximation is that Z,,(O) = 1, hence the first E.-moment has an asymptotic behavior as l/r at small r, matching that of the scattering irradiance. As a consequence, this implies that (4,)
130
JOHN LANGERHOLC
will not approach 0 as r + 0, but rather lim (4,) = L r+O
(8.25)
cwJ?r.
Calculation of the Second Moments For the second moments, the quantities to be inserted into the inverse transform are given by 8.6.
Qz,o(l))=f2,o(ll)j(q,~fl~=0o),
(8.26)
Qo.,(~)=fo,~(‘l)j(tl,zfl~=~),
(8.27)
where
- (?(I-exp(
-(nr,i20)2))iL. (8.28)
fo,e(
11) =
-
f,(9
2P2
-(
q2-
7x2
1
=-+
tly
-“(O, 2 ffsqs
F(l-exp(
f,,
y(o,
7,
q)
8, x) + $$exp(
- (qz//2ar)‘]
-(qz~/2*)2]])2.
These each decompose into three integrals of different character. The constant terms l/2/3’ merely reproduce a multiple of the irradiance and need not be calculated separately. The remaining two are of the types [ Q.21 and [ Q.31 ,involving 119and T$, respectively. The formulas for the angular spread may then be written
(Ab)2=l/~z+{M,,(r)+ M&-)}/i(r-,q)
(8.29)
131
Forward Scattering Amplitds where
i%-,(v) =
2 z ( )I
ufzf
(y2
ew(- (vq/2a)2)
- y{l-exp(
-(qzf/2u)2])2
(8.30)
f(r,zf) I
(8.31)
(8.32)
2 OfZf
( ii
exp( - (wf/2~)2}
-f-$(1-exp{
-(?JZf/201)2]]2
A&,(rl)= :
(112
i(r,zf).
(8.33)
I
The integrals with
involving
7: are most conveniently evaluated using [Q.3]
(8.34)
rlzf ,2a)2])2]
i(r,
zf),
(8.35)
132
JOHN LANGERHOLC Fexp(
9ee( V/T) = c
-(tlzf/2cYq2)
i -~{~-~~~{-(n~f,2~~)2)}2
eucs),
(8.36)
I (8.37)
we find
A4&-)=(C/2*r%2 'J - 2u$zr2
0
O1 afZfexp (-(l)zf/2Ur)2]
(8.38)
i
[I-exp(
-(qzf/2ar)2]]2
e”(@J,(v)dv. I
In the second integral the factors of J1 fall off either exponentially fast as tm2. To evaluate MO,, we apply [Q-2] with n = 2 and
9&)
= - ‘(z2zf)
i( r, zf)
or as
(8.39)
and use the fact that dJ,(q)/dq = - J1(q) to obtain
M&r)
= (C/2rr2az
)[av~(n)eu(v)41(n).
(8.40)
Sufficient falloff can be obtained by integrating once by parts
Mar(r)
= - (C/2rr2a2
) L=Jb) dq
~{~(1))+r)~‘(rl)+rl~‘(rl)~(rl)}e~(”). (8.41)
Forward Scattering Amplitudes
133
To use [Q.2] to evaluate M,,,
4rr(9) = -
for which n = 2 and
exp{ - (rIz//2o)‘)
y
1 - 7
{ I-exp(
- (qzf/20)z]]2)
i(r,
zf),
(8.42)
we rewrite, from Bessel’s equation (8.43)
-~on(rl)=Jo(rl)+J~(rl)lrl=lo(9)-J~(rl)lrl. Insertion yields, finally,
M,,(r)
= ( C/27rr2a2
x
urzfexp [ - (11+q2] i
- 2u-J2r2
(1-exp{
eWcrl). (8.44)
-(~zf/2ar)2])2
1 9.
RESULTS
The Gaussian parameters were computed as functions of the radial distance r from the beam axis in the focal plane zf = 2 km for a laser of F = 1 megawatt for several values of scattering parameters in the wasp-waist approximation. The maximal irradiance values were estimated by the Gaussian fit and plotted in Figure 3 for the values ur = 0.035 km-’ and o = 15 radd’. Owing to lack of data on the absorption coefficient, the extinction was arbitrarily set at ue = 2 a, in all computations. The values can be adjusted for different values of the absorption coefficient by multiplying with the appropriate exponential. It is seen that the maximal radiance value is exceeded within 4 meters of the beam axis. Figure 4 shows a plot of the single-ray Green’s function for the irradiante for the values uf = 0.035 km-’ and (Y= 8.11 rad-‘. Figure 5 shows the same function for the larger value of uf = O.l6/km. In spite of the
JOHN LANGERHOLC
134
(W/m%) 4
lOe?--
lOe6-R d" ' a n c e
lOe5--
10a4--
10e3
! 0
b
I
2
4
6
6
10
12
14
16
16
20
Radial Distancefmm Axis
FIG. 3.
(Wlm2)
I
0
50
l (m) 100
150
200
250
Radial Distance from Axis
FIG. 4.
300
350
Cm)
135
Forum-d Scattering Amplitudes (Wlmz) 1000.00
1 I
r
r
100,oo --
d”
I
a n
IO,00 --
C e
I,00 3
0
5
10
15
20
25
30
35
40
45
50
ä(m)
Radial Distance from Axis FIG. 5.
increased absorption coefficient the values of Z are approximately three times as great. The angular spreads A+ of the distributions are generally small. They are less than 0.1 rad at r = 10 m for both values of the scattering coefficient u and in both r and 8 directions. performed numerically in the The power integr ai s 2 rr / Z(r, zf)rdr focal plane up to r = 370 m agree with the theoretical values integrated over the entire focal plane to within 4% in both cases. The author is indebted to G. Born and H. Hermansdorfer of MesserSchmitt-Bijlkow-Blohm, 8012 Ottobrunn, West Germany, for financial support of this work, and especially to G. Born for suggesting the problem and its relevance to questions of laser safety. He wishes to thank M. Gorriz, also of MBB, for discussions and assistance with the computations, and P. Glijckl for the figures.
APPENDIX.
LIST
OF
SYMBOLS
Z(+,r, z): The radiance
as function
of the angular orientation
coordinates (r, 2). Z(r, z): The irradiance, the spatial coordinates
or angular (r, z).
integral
of the radiance,
4 and spatial as function
of
JOHN LANGERHOLC
136 f(f, 1, z):
The
Fourier
transform
of the
radiance
in angular
and radial
variables.
f(q, z): The Fourier = i(O,q, 2). CY:Inverse breadth
transform
of the irradiance
of the forward
in the radial
variables
peak of the phase function.
ur: The coefficient of forward scattering. (Y,: Inverse breadth of the reverse peak of the phase function. a,.: The coefficient of backscattering. a,,: The coefficient of isotropic Us: The absorption coefficient. a,: The extinction coefficient. p: Inverse y: Inverse
scattering.
angular divergence of the input beam. beam radius RfoL in the focal plane z = zf.
Rfoc: Beam radius (l/e) in the focal plane x = zf. Ri,: Beam radius (l/e) in the input plane. zf: Distance from the input plane to the focal plane.
REFERENCES 1 2 3
4 5 6 7 8
9
10 11
H. S. Snyder and W. T. Scott, Phys. Reo. 76:220-225 (1949). L. R. Dolin, Propagation of a narrow beam of light in a medium with strong anisotropic scattering, Zzv. Vyssh. Uchebn. Zaved. Radiofiz. 9:61- (1966). R. L. Fante, Propagation of electromagnetic waves through turbulent plasma using transport theory, ZEEE Trans. Antennas Propagat. AP-21:750-755 (1973). A. Ishimaru, Wave Propagation and Scattering in Random Media, Academic Press, New York, San Francisco, London, 1978. A. Zardecki and W. G. Tam, Laser beam propagation in particulate media, Can. 1. Physics 57:1301- (1979). W. G. Tam and A. Zardecki, /. Opt. Sot. Am. 69:68-70 (1979). W. G. Tam, Multiple scattering corrections for atmospheric aerosol extinction measurements, Appl. Optics 19:2090-2092 (1980). W. G. Tam and A. Zardecki, Multiple scattering corrections to the Beer-Lambert law. 1: open detector Applied Optics 21:2405-2412 (1982); 2: detector with a variable field of view, A&. Optics 21:2413-2420 (1982). A. Deepak, V. 0. Farrukh, and A. Zardecki, Significance of higher-order multiple scattering for laser beam propagation through hazes, fogs, and clouds, Applied Optics 21:439-447 (1982). J. Langerholc, A beam-width formula derived from a two-component scattering theory, in preparation. J. Langerholc, Volumes of diced hyperspheres: resumming the Tam-Zardecki formula, Appl. Math. Comput. 30:1-18 (1989).
Forward 12
13
14
Scattering
Amplitudes
137
‘K. Altmann, Forward-scattering formula of Tam and Zardecki evaluated by use of cubic sections of spherical hypersurfaces, Appl. Optics 28:4077-4087 (1989). J. A. Weinman, J. T. Twitty, S. R. Browning, and B. M. Herman, Derivation of phase functions from multiply scattered sunlight transmitted through a hazy atmosphere, J. Atmos. Sci. 32:577- (1975). A. Zardecki and S. A. Gerstl, Multi-Gaussian phase function model for off-axis laser beam scattering, Appl. Optics 26:3000-3004 (1987).
‘Although perhaps now only of academic interest, this paper could have been considerably shortened and the unnecessary approximations avoided by using the original recursion relation of [ll] to compute the volume functions iteratively. In fact, the area functions themselves need never be computed if a Stieltjes integral is used, since the appropriate measure Surf.(r) dr = dVol.( r) is in fact the differential of the volume function.
Strassbergerstr. 10 8000 Miinchen 40 Munich, Germany
ABSTRACT The equation for radiative transfer is solved in the case of forward scattering of a laser beam that may be focused into a medium to attain its minimum in uacuo width beneath the surface. The moments of up to second order in the angular distribution of the radiance are calculated for an arbitrary point within the scattering medium and used to obtain a Gaussian fit to the radiance distribution at that point in order to reduce the mass of required computation to a small set of fundamental parameters essential for considerations of design and safety. The Gaussian parameters -u-radiance, tilt, and angular breadths of the radiance distribution-are Fourier transformed back to configuration space, where they appear as one-dimensional Hankel transforms. For calculations of the off-beam scattering amplitudes in the focal plane of the unscattered beam, a “wasp-waist” approximation is introduced, in which the radial profile of the beam shrinks to a S-function in the focal plane. A method for the improved convergence of the Hankel transforms used to evaluate the Gaussian parameters in this approximation is derived. It is then shown that these single-ray solutions, which are mainly independent of the details of the (narrow-angle) beam, can be used as Green’s functions, enabling the exact scattering amplitudes to be recovered by convolving them with the normalized profile of the unscattered beam in any plane of interest. Thus even for broad beam profiles and points near the beam axis, the computation of the exact Gaussian parameters requires only a single further integration of the relevant Green’s functions over everywhere positive, i.e., nonoscillatory, functions. The asymptotic behavior of the Green’s function for the irradiance is studied at small radial distances r, and an approximative formula for this function is derived.
1.
INTRODUCTION
It is well known that atmospheric scattering degrades the performance of laser systems, reducing visibility, attenuating the direct rays, and blurAZ’Z’LiED MATHEMATICS AM) COMPUTATZON 48:101-137 0 Elsevier Science Publishing Co., Inc., 1992 655 Avenue of the Americas, New York, NY 10010
101
(1992)
0096-3003/92/$05.00
102
JOHN LANGERHOLC
ring contours. In atmospheric applications of high-energy laser beams, scattering can not only reduce the amount of power delivered to its goal, but can even pose a safety hazard by rerouting significant amounts of energy outside the intended path of the beam. This is true not only for observers situated near the axis in the path of the broadened beam and looking in the direction of the transmitter, but also for observers at the transmitter site looking down the length of the beam. The relevant effects of back scattering from a focused beam will be studied in a follow-up work; the foundations are laid here in the calculation of forward scattering amplitudes at arbitrary depths. Although the solution of the transport equation in the small-angle approximation has been known for some time [1-lO],previous investigators have experienced numerous difficulties in numerical evaluation of the inverse Fourier/Hankel transforms involved in off-axis calculations of forward scattering amplitudes. Various methods have been suggested for evaluating the fourfold inverse Fourier transform required to determine the angular distribution Z(&r, z) for a given point in space (r, z). A Taylor series expansion of the radiance Z(&‘,r, z) about c$’ = 4 in the convolution term has been proposed by Fante [3]. Significantly different results were obtained by Tam and Zardecki [6] by expanding the exponential en in a power series and integrating term by term. Even in the relatively simple case of a conical beam diverging with Gaussian angular spread from zero width at the origin, however, the computational effort required for the evaluation of the nfold integral occurring in the nth term proves to increase so rapidly that the latter authors break off the summation when the last summand becomes less than 5% of the total attained. This results in taking from 3 to 6 terms, depending on the axial distance. The resultant expressions for the nth term are sufficiently complicated even in the calculation of the irradiance, which no longer depends on the angle 4, that special methods have been devised to improve convergence [ll, 121. In [II] ,a resummation of the series was achieved, leading to an integral over a K, function with exponential falloff, but complex argument. In [12],summation was carried out by a rather complicated computer program. In the light of the present, much simpler results, neither of these procedures appears very satisfactory. 1.1. OjfBeam Zrradiance In addition to the difficulties already experienced for a conically diverging beam, the desire to calculate in the waist plane of a convergent, sharply focused exacerbates the problem touched on in [9] ,where
in such calculations off-beam amplitudes laser beam further the kernel r~Ja(r~r)
Forward Scattering Amplitudes
103
oscillates with increasing amplitude and is only very slowly damped in the integration variable 7 for radial distances r larger than the beam width. 1.2.
Data Compression
To calculate at a given point the radiance, which depends on both the position of the point and orientation, it is theoretically necessary to evaluate a four-dimensional inverse Fourier transform. For many considerations, including those of laser safety, where it is a question of maximal exposure, it is not necessary to calculate the exact functional dependence of the radiance on the orientation angle in those cases where the breadth of this distribution is so small as to render all the light scattered to a given point effective. This condition is expected to be met in the case of narrow-angle forward scattering, and in the case of sharply peaked back scattering to be treated in a subsequent work. For this reason, a theoretical approach will be followed here which first reduces the computational volume from a fourfold inverse Fourier transform at each point (r, z) and every orientation 4, to a set of twodimensional transforms, condensing the essential information into a set of relatively few parameters basic to the angular distribution, including irradiance, average orientation and angular breadths in the radial and tangential directions. Since these measures figure as parameters in a Gaussian representation of a beam, they are referred to here as the Gaussian parameters of the scattered radiance distribution. In cases where these parameters in themselves are not sufficient, the radiance can still be estimated as a function of the orientation angle Q, by calculation from the Gaussian distribution determined by the calculated parameters, as discussed in the following subsection. 1.3. Focused Beam The configuration to be studied is illustrated in Figure 1, showing a cylindrically symmetric laser beam focused into a scattering medium, where it attains a minimum width of Ri, = l/y at a penetration depth z = zr and then diverges again at greater depths. Although originally intended for this special case, the calculations are general in the sense that zr can be chosen negative or zero to place the focus outside the medium or at its boundary. The latter case is the one that has been studied almost exclusively up to now. The scattering phase function will be assumed as in [5-g] to have a Gaussian forward peak; residual isotropic scattering will be treated in subsequent works by means of a two-component scattering theory [lo].
JOHN LANGERHOLC
104
-----
_-._,-.- .._.-.- .--. -.-.
Input
FOCd PlalW
Plane LO
z=z, FIG. 1.
2.
RADIATIVE
TRANSFER
THEORY
The basis of the present calculations will be the radiative transfer or transport equation derived according to geometrical optics from energetic considerations. Interference and polarization effects are neglected; the fundamental quantity to be solved for is the radiance or specific intensity distribution Z(+,x). The solution of the basic equation in the small angle approximation will be reviewed to assist in generalizing the known solutions to the case of a beam not necessarily focused at the input plane. The Transfer Equation The basic equation for the radiance derived from consideration of a differential cylindrical volume element with generator in the direction of the unit vector 1, is [l-4] 2.1.
l,*v,q~,x)+ uJ(4,x) = Uf// +wy(4’,x)
dfl,,.
(2.1)
Here P(+, 4’) is the ph ase function, giving the probability that a ray will be scattered from the angular direction 4’ to q5. It is normalized to unity: s 4a
P(&,&)
dn,
= 1;
(2.2)
and the scattering coefficient is designated by ur, the amount of power diverted out of a differential volume dV = dA* dz being given by = ufzdv= dZSCatt
(ufdz).(Z-dA).
(2.3)
105
Forward Scattering Amplitudes
Similarly, the extinction coefficient gives the total power lost from the beam when traversing the elementary volume: dP’=
a,ldV.
(2.4
In the multicomponent scattering theory of [lO],it is composed of the various scattering coefficients corresponding to forward peak, isotropic component, and rear peak, along with the absorption coefficient: u,=uf+a,+u,+u,.
(2.5)
For the present work dealing with forward scattering, two terms in this equation will be relevant.
only the leftmost
2.2. The Small Angle Approximation General solutions for the radiative transfer equation are difficult to find; and it is customary in the case of forward scattering of a collimated laser beam to make the following approximations, based on the assumption that the bulk of the scattered light propagates in directions deviating by only a small angle from the principal axis 1 z. of the laser beam [l-9]. The coordinates are represented as x = (r, a);
and the angle C#is spanned by a two-dimensional
(2.6) vector
1,=(&l),
P-7)
where
The angular integration thus turns into a planar integral:
In the remainder of this paper, the limits of integrals that are taken from - 00 to 00 will be omitted. Integrals without limits should be understood as such.
106
JOHN LANGERHOLC
The phase function is also considered to depend only on the magnitude of the vectorial difference of incident and scattering angles:
p(w)= With these approximations becomes
=U
f
q M-4’1).
(2.9)
and changes of notation, the transfer equation
JJ
P(+ - ti’)I(+‘,r,
(2.10)
2) d24’.
Fowi.er Transformation of the Transfer Equation 2.3. In this form, the natural method of solving the equation is to perform a fourfold Fourier transform in the coordinates 4 and r leaving the axial coordinate z as an independent variable:
~^(t,q, z) = J
J
J
J
(2.11)
d2~d211e’(0’~+II”)Z(~,r,z).
The difficulty of solving the equation in physical space is thereby transformed into the difficulty of evaluating the fourfold inverse Fourier transform to obtain the radiance distribution:
z(4,rF 2)=(27r-4 J J J J
d2~d2rle-‘(E-~+?‘r)i(f,tl,
z).
(2.12)
It turns out not to be advantageous to transform the equation in all variables at once, since the factor C$ would turn into the differential operator - i V{. 2.4. Fourier Transformation of Geometrical @tics The basic idea of the solution is to transform the ray propagation condition of free space Ar = +Az into a substitution for the conjugate variables which will allow the solution to be extended from the input plane z = 0 by a simple variable substitution. The equivalent of the ray propagation condition will turn out to be the substitution
E+f+zv
(2.13)
107
Forward Scattering Amp&&s
inserted into the transformed radiance at z = 0. This will not only allow the solution to be written down relatively simply from the transformed boundary values but will also provide insight into the problem of representing, and solving the transport equation for, a focused beam. 2.5. Spatial Fourier Transformation To put this statement on a rigorous basis, the transport equation is transformed at first in the two radial coordinates only, giving [4]
(2.14) The constant factors of the unknown on the left-hand side can be converted into an exponential factor by introducing the intermediate unknown I’, which is defined by:
i(4,v,r)=exp(- o,z+i+-rlz)i’(+,q,z)
(2.15)
This simplifies the equation to
13,i’(c$,q, 2) =
uf//
P(+- q5’)e-i(“-“‘)‘qZ
x i’($‘,q, z) &jJJ. 2.6.
Solution of the Transport
(2.16)
Equation
Performing now the Fourier transformation in the angular variable and using the convolution theorem, we find
where (2.18) The solution of the last differential equation can be immediately written down; it is:
Q,% z) = e”(f*-~*t)?(f,tj,O),
(2.19)
108
JOHN LANGERHOLC
where
and
qt,rl>q = @,%O). From
the above definition
of I’, it follows that
=exp(-o,z)i’([+qz,q,z). Inserting
the solution obtained
(2.20)
for Z’ yields
i(l,g,z)=exp(-~~z+62(S+11z,-rl,z))i’(~+?z,r1,0). Using the equality
(2.21)
of Z and I’ for z = 0 and the fact that
n(o,?z,-~r,~)=u~J’~(E+a(z-z’))dz’ 0
=(J
f
J
‘p(E+Ilz’)dz’=n(E,~,z),
(2.22)
0
we have, finally, i($,~.z)=exp(-
a,~ + Q(E>l),
z))Z(f
+ MLO).
(2.23)
This may be expressed in words as follows. On propagating through a depth increment z, the Fourier transform of a given input radiance distribution is changed in three ways. The conjugate variable E of C$ is replaced by t: + gz; and the entire expression is multiplied by the extinction factor and by e”.
Forward Scattering Amplitudes
109
Contained within this general prescription for translating the Fourier transformed radiance in the z direction is a substantiation of the claim that the radiance of a beam propagating in uacw, can be calculated from its Fourier transform at z = 0 by replacing the conjugate angular variable [ by E + Qz. It will be used in the following discussion to derive a representation for a beam focused into a medium. This behavior is also reflected in the argument E + 7 z’ occurring in the Fourier transform of the phase function. 3. 3.1.
THE
PHASE
Representation
FUNCTION of the Forward
Peak
The conventional representation of a forward peak in the phase function is a normalized Gaussian in the transverse components of the (tangent 04 the angle [S-9]:
P( 4) = :exp(
- (~‘4’)
JJ
P(+)d2$=1.
(34
The Fourier transform is also Gaussian: F(E)
=exp( - t2/4a2).
(3.2)
This representation of P leads to the following expression for the exponent of the scattering function:
Q(e,v,~)=u~~‘exp(which is expressible
(f+qz’(‘/4cY2)dz’,
(3.3)
in terms of standard functions:
X (erf(
nz+~~se)
-ej
F)),
(3.4)
110
JOHN LANGERHOLC
where 0 is the angle between the vectors t and q:
P-5)
t’q=[qcose. For [ = 0, this simplifies to
(3.6) 3.2. Multi-Gauss Peak In the case of certain phase functions, it is desirable to take a multiple forward peak into account [IS-141 ,in which a multi-Gaussian phase function is constructed as a convex linear combination of single Gaussian peaks of different widths: P(f#l) = r-l
5
exp (-
K,c$
m=l
c&“)>
(3.7)
where
5
K,,,=l.
m=l
The result on the solution for the scattering amplitude is that Q is replaced by the same linear combination of the Q functions calculated for the individual scattering widths:
P-8) 4.
REPRESENTATION
OF THE
INPUT
BEAM
To find the forward scattered radiance from a laser beam directed into the medium, it is necessary to insert the Fourier transformed radiance of the beam at the input plane z = 0. Gaussian Beam The radiance distribution for the laser beam used in [9] is taken to be Gaussian in both angular and radial variables in the input plane
4.1.
Z(f$,r,O)
= FP2y2ap2exp(
-
P2+2- y2r2),
(4.1)
111
Forward Scattering Amplitudes
where F is the total input power of the beam. The Fourier transform is again Gaussian: i([,q,O)
= Fexp(
-
52/(2fl)2-71~/(2’~)‘]-
(4.2)
According to the general solution found above for the radiative transfer equation, the solution propagates in vacua according to the substitution
i.e., in the absence of scattering and absorption, the radiance in the plane at depth z is given by I^(~,ll,~)=Fexp(-(f+112)~/(2p)~-8~/(27)~).
(4.4)
4.2. Inward-Focused Beam The input beam described in the preceding section has its minimum breadth W(0) = v&/y in the input plane z = 0 [3, lo]; it spreads with increasing z, having a width at any given depth of
R(z)”= iw”( z) = l/y2
+ ( z/q2.
(4.5)
This solution is symmetrical about the input plane z = 0; it represents a beam converging to a relative focus at z = 0 and then diverging at the same angle l//3. To represent an unattenuated beam focused at a given depth zr beneath the surface, it is merely necessary to translate this unscattered beam along its own axis by substituting z - zf for z:
P-(&q,
z) =
Fexp{ - (t + II( 2 - ~f))~/(2/3)“gels)-
(4.6)
Thus the input radiance of the infocused beam is
Pm(t&o)
= Fexp{
- (E - Irzf)“l(2b)“-
q2/(2y)“].
(4.7)
Its width at the input plane is given by
R(0)”= iw”(o) = l/y2 + ( zf/q2.
(4.8)
JOHN LANGERHOLC
112
By means of these relations the distribution can be fully determined by the widths in the focal and input planes and the separation zf between the two planes: r=JUw(z,),
p=
v%Zr/J{W2(0)-w”(q)} . from a
Radiance Scattered
4.3.
(4.9)
We
are now in a position Fourier space for the scattered the input
plane
as boundary
(4.10)
Focused Beam
to write down the radiance by inserting value
in the general
complete solution the beam radiance solution
found
in at
above.
The result is
f(l,q, z) = Fexp(Xexp{ 4.4.
Q(~,T, z))
uez + -(I
+ v( z - 9))“/(2S)‘-
v~/(~Y)‘}.
(5.11)
Representation of the Beam in Configuration Space
The inverse Fourier transform of the unattenuated beam at a given depth z, measured again from the focal plane, can be performed analytically after the variable
substitution:
d2T = d2f
7=c+7jz; which
acts in the exponential zbeam(‘pr,
z)
=
to separate
the integration
(2T)-4F/- 1 dz,,/ J X exp( - r2/4P2
= Fb2y2T-2
(4.12) variables:
d27e-‘(‘.*+‘.-““.&)
- q2/4y2)
exp( - (/3+)‘-
With r and z constant, the &dependence Gaussian form by completing the square:
y2(r-
2#)‘}.
may be explicitly
(4.13) rewritten
in
~82+(~y)2)~2-2zy2r~d=(P”+(~y)2){~~y2rl(B2+(w)2]]2
- y4z2r2/(
P2 + ( z~)~].
(4.14)
Forward Scattering Amphuh Insertion
113
above leads to the expression
where
Zy)2)= l/~2’Y2R2( z,
(A+)” = l/( f12+(
zry2=
zry2(A+)‘,
(4,) =
P2 + ( zY)2
R”( 2) ‘iW2(
z) = ( z//3)2+l/y2.
This may be taken
as a motivation
dependence
radiance
of the
for the Gaussian
distribution
fits to the angular
to be developed
in the
next
section.
5.
GAUSSIAN Rather
than
FIT
TO THE
producing
ANGULAR
tables
DISTRIBUTION
of numerical
values
by means
of an
intensive computational procedure that evaluates the radiance in every angular direction + for a given point (r, z), the orientation, magnitude, and widths of the bell-shaped angular distribution will be calculated given point. This will allow the functional form to be modeled equivalent two-dimensional Gaussian distribution of the form
=
-
W)A4
?rA+
r
.9
exp( - (4,
-(+,))“/A@
- +82/A&),
at the as an
(5.1)
and enable an adequate fit to the distribution. The physical significance of the parameters is as follows: Z(r, z) is the irradiance, the angular integral of the radiance at the given point. The mean or expectation value ($,) of the radial component of the orientation angle 4,. is a measure of the tilt of the radiance lobe outward along 1, away from the principal direction of propagation of the laser beam 1 t. The
114
JOHN LANGERHOLC
corresponding peripheral component of the tilt vanishes by symmetry. The two widths of the distribution A+, and A+, are analogous to the inverses of the parameters CY,P, and y for the phase function and the input beam. The advantage of such a Gaussian representation is that the parameters are independent of the orientation angle 4; hence only two-dimensional Fourier transforms need to be evaluated, rather than the numerically more demanding four-dimensional transforms of the radiance function. 5.1. Moments of the Angular Distribution The fundamental parameters of the fit can be calculated from the moments M m,“[ I] up to second order (m + n < 2) in the angular variables, where Mm,,[f]
-/
/
(~~)“(~y)“f(~)d2~.
(5.2)
To avoid unnecessary notational complications, we take advantage of the cylindrical symmetry of the problem and align r along the r-axis by setting y = 0. In so doing, 4, becomes 4, and & becomes c$,,. The zeroth order moment is simply the irradiance:
M,,,[ I] = /
1 Z(hr,
It is obtained by inverse Fourier
Z(r, z) = (27rP2/
z) d24 = Z(r, 2).
(5.3)
transformation / e-““‘Z(q,
z) d2q
(5.4)
of the solution obtained above
i(% 2) = i(0,rl.z). The expectation values of the angles are obtained corresponding first moments by the zeroth moment:
(4%) = M,,o[11/%lJ[ 11. As remarked previously, (&J will vanish by symmetry.
(5.5) by dividing the
(5.6)
115
Forward Scattering Amplitudes
If the Gaussian fit I, is inserted for I in the calculation of the second moment, we have
x exp( - (4, -
(~J,))~/W~ - 4,‘/&,‘] I(r, 2)
= { (4,)”++br)‘}Mo,o[~] . Thus the condition for matching the second radial moment, M,.,t
(5.7) M,,,[ &] =
Zl,yields +4r)”
=
M,,o[
11
Po,o[
I]- (4,)“.
(5.8)
Similarly, but more simply,
,(A40)2= ~0,2[~1/~0,0[~1-
(5.9)
5.2. Moments in the Fourier Representation If the Fourier transform of the distribution is known analytically, the moments are readily calculated by substituting - iV[ for 4 and evaluating the appropriate derivative of the distribution at f = 0, which corresponds to integrating over all values of 4. Hence the moments of interest, still Fourier transformed in the spatial coordiante r, will be given by M,,,[ i] = - ial(fJ,tl, M,,,[
~)/a,$,,
i] = - d21(0, 1, ~)/a[,“>
M,,,[ i] = - a21(O>~, ~)/a$
(5.10) (5.11) (5.12)
To obtain these quantities as functions of the spatial coordinates, it suffices to perform the two-dimensional inverse transformation in 7. 5.3. Di&mntiating an Exponential The relevant part of i for the differentiation in E can be represented as ef(z), where f(E)S-{t+7(z-
2/)}2/482+n(t,a,z)+Const(r,z),
(5.13)
116
JOHN LANGERHOLC
and where the last term does not depend on 1. Because of the exponential form, we may write aefcQ/i3E, =
fx(t)ef([),
(5.14)
t )) efl%
Peflf)/8Ei
= { fX( f)2+
f,, J
Cr2ef(l)/a[i
= ( fY( e)‘+
f,, &f
))@).
(5.15) (5.16)
To obtain the moments as functions of the spatial point (r, z), it will then suffke to insert the appropriate differential expressions into the Fourier inverse transform integral for the irradiance, the correspondence being:
$-$I, 1, 2)
M 1.0
*
M 2.0
--f*(o,l,2)2-f,,.(9rl,z)
M 0.2
--f,(o,r,,r)“-~~,,(3,,,r)
-
(5.17)
With the expression given above for f, we find f,(e)=+,++-
~~)}/2P~+an(~,rl,z)/aE.,
(5.16)
where
At f = 0 the last integral can be evaluated analytically:
aqo,tl, 2) = - 7 at,
J,‘exp{
-(qz’/2o)“)
d(gz’/2a)2
(5.19)
117
Forward Scattering Amplitudes
For the first derivative, then, we have ,(,)JJz-zf) x
-~[l-e”p(-(~z/201)2)).
2P2
(5.20)
The second derivative is f,,,(E)
= -1/2P2-
Q(k%
z)/202
+ qzz’)2exp[
- (4 + ~z’)~/~cY~} dz’.
(5.21)
The last term can be simplified for ,$ = 0 by partial integration of
J‘(
qx
z’)‘exp(
-
( 1Z’/2~)2]
dz’
0
= -
exp( - (qz’/2o)‘)
di-
zexp[ - (qz/2a)2]
so that
- gexp{
(5.23)
- (qz/2cu)2).
Similar results are obtained for the derivatives with respect interchanging qX and vy. 6.
GAUSSIAN FIT PARAMETERS
AS HANKEL
to t,
by
TRANSFORMS
The quantities to be calculated for a given vectorial separation r from the beam axis at the penetration depth z have been derived above as
118
JOHN LANGERHOLC
functions of the Fourier following form:
conjugate
variable
q, and can be written in the
where i is the Fourier transform of the irradiance i(‘l.z)=F.expI-Sr1’/R2(z)-u~~+n(?,~)},
(6.2)
and R(z) is the beam width at depth z: R”( 2) = (2 - z_J2/4p2 + l/y?
(6.3)
The transformed n-radiance depends only on the magnitude of q, whereas the f(q, z) may contain combinations of r,rr, 7: and $. The asymmetry arises because r = (x, 0) is chosen parallel to the r-axis so that r = X. The quantities Q, including irradiance, tilt, and angular breadths of the radiance distribution, are calculated as inverse Fourier transforms of the Q:
Q(r? 2) = (2~)-“/
/
d”?e-““.‘@~,
z),
(6.4)
where
d2q = qdvde. In the sequel, for the sake of simplicity, the calculations will be performed in the focal plane z = zf. Hence the depth coordinate will be understood to be zf and suppressed in the functional arguments of Q. Cylindrically Symmetrical f If Q is independent of 8, the angular integral becomes a Bessel function of zeroth order: 6.1.
I
0
2*
de e-iqr.cose = 2rJo(53-),
(6.5)
119
Forward Scattering Amplitudes
hence
Angular Dependent f If &II)= (- irlx /~)“drl), then
6.2.
JJd211e-i1).'6(1)=(d/dx) J
= p/w
J
d2ve_iq’rg-“q(q)
J J 7l~a9(~)a-“Jo(v)
and hence
If, fin&,
rewritten
J
= 9(v) sin28, then the O-integral can be
O(v) = (71,/9)~9(9)
2* e
-'*"OSe,in2ede=(i*r)lJo26ine
&-ivcos0
0
= (i/qr)
de=2?rJ,(vr)/qr, J2=e-isrcoseCOSe
(6.6)
0
and hence Q(r)
J J d2~e-iq.rO(~)
= (2r)-"
=-
1
2*?-
=& Joo~~~++7(+-) [QJI Jmds./l(sr)s(v) 0
120
JOHN LANGERHOLC
6.3. Fallofl Properties of the Q(v) Difficulties in calculating the Jo-transform for the irradiance were encountered by previous researchers [5-91 because of the extremely slow falloff of the function &q/r) = i(v/r, z) in the variable 7. Combined with the factor v arising from the measure d2~ = r] dv de and the only very slightly damped oscillations of the Bessel function, this means that a direct evaluation of the integral
P-7) as it stands will involve integration over a large number of lobes of alternating sign, which initially increase rapidly in absolute value (as Jq). The method is recommended in [9] only for radial distances from the beam axis on the order of the radius of the unscattered beam or less, since the factor exp( - +q2R2( z)) will provide sufficient falloff in these cases to limit the required range of integration, hence the growth of qJa(s) within the range of interest. 6.4. W-Beam Zrradiance For such applications as calculations in the focal plane, this state of affairs is quite problematical, however, since the irradiance and other parameters of the Gaussian fit are to be evaluated at radial distances for which the unscattered beam has practically zero amplitude. After substitution of q/r for q, the falloff contribution from the beam width becomes which is practically ineffective for r] of moderate exp( - $R2( z)/4r2), size. The situation becomes dramatic if at the desired depth zf, the focused beam attains a minimum unscattered radius R( z,-) = l/y on the order of 1 cm, for instance, whereas a typical radial distance of interest would be on the order of r = 10 m. Falloff from the damping factor begins to become appreciable only above q = 2 r/R(O) = 2000. 6.5. Excision of the Unscattered Beam Since its value is not only known, but in this special case is expected to be practically zero, the contribution from the direct beam need not be calculated along with the scattering amplitude, but may be subtracted from the Fourier transform of the irradiance before evaluating without materially changing the value of the result to be computed. This value is obtained by setting Us = 0, hence Q = 0, in the above expression for i, while leaving the original extinction coeffkient a, unchanged.
121
Forward Scattering Amplitudes In the remainder
of this calculation,
tion to the irradiance
thus, only the scattering
contribu-
transform
il(~,z)=F.exp(-~R2(a)r12-u~z)(e”(?.Z)-1)
(6.8)
will be transformed back to configuration space. Since qQ(q, z)+ aa,& as ~+oo, the last factor on the right hand side will fall off as l/q. After the substitution
of q/r
for 7, the exponent
becomes
n(q/r, 2) = (ckwfrJ;;/7j)erf(7pz/2ar). The falloff due to the factor (e” - 1) can be estimated, following inequality: lea-11
(6.9) if desired,
from the
***)
(6.10)
The Wasp- Waist Approximation
6.6.
In view of the numerical facts just discussed, the factor exp( - aR”( zr)v2), which provided the sole damping in the original expression, can now be completely ignored, since within the effective range of the Bessel function, which falls off as l/~/q, its value will be essentially equal to 1. The convergence of the integral is insured by the fact that the maximal values of lo oscillate and decrease monotonely in value to zero. The result after subtraction
is invariant with respect
to the details of the
beam $(r]/r, as could be expected
Z) 2: F*exp(
- u,z)(e”(““XZ)-l),
from the small angle approximation.
(6.11) These
details
will not figure in the angular integrated irradiance but, by virtue of their appearance in the f(q/ r ), will enter into the results for the tilt and angular breadth of the radiance distribution. Geometrically, this approximation is equivalent to assuming a conical beam with radius R( zr) = 0 in the focal plane. For a sufficient radial distance within this plane, the constriction of the beam will have practically no effect on the amplitude of forward scattering at distances large with respect to its focal width. (The calculated value at r = 0 will nevertheless diverge, clearly showing the limit of this approximation.)
122
7.
JOHN LANGERHOLC
RECOVERY
OF THE
EXACT
IRRADIANCE
The usefulness of these calculations is not, however, limited to the focal plane or even to calculations of an inward focused beam. On closer inspection, it can be seen that it need not even remain at the level of an approximation. In the sense of function theory, the omitted factor exp( - aR”( 2)~~) represents a convolution of the inverse transformation of the single ray factor (en - 1) with the normalized unscattered beam profile in the plane at depth z. For off-beam values of r, where this function varies slowly with respect to the radius of the beam profile, it represents only a negligible averaging (or smearing in the sense of distributions) of the Green’s function from the single-ray solution and thus has a negligible effect. For calculations near the axis or at penetration depths for which the beam radius is appreciable, the error made in the approximation can be undone by subsequent convolution of the single-ray solution lr in configuration space with the reduced beam intensity normalized to unit power flow:
II(”
JII(r’.ZIR=O)IG(r-r’,z)dZr’;
(7.1)
where
&(r, z) = Z(r, z 1uf = O)/F( 2) = exp( - r2/R2( Z))/7rR2( z),
IS
(7.2)
Zb(r, 2) d2r = 1.
This procedure is illustrated in Figure 2. Here again it is understood that to obtain the reduced beam intensity I,, the extinction coefficient ue is maintained at its original value; the scattering coefficient ur is set = 0 as an independent variable. In the case of a Gaussian beam profile, the angular integration can be performed analytically. It leads to a convolution of a product of a Gaussian
123
Forward Scattering Amplitudes
FIG. 2.
with the Bessel function of imaginary argument I,:
(‘-3) where
i( r, I-‘) = f
s
* dtl Zb(r-r’, --*
z)
dBexp{2rr’cosO/R2(z)}
(7.4)
124
JOHN LANGERHOLC
Once the Green’s function Zf(r, z 1 R = 0) is calculated for all values of r, the total scattering amplitude for any axial distance from the beam axis is obtained by a single one-dimensional integral over a positive integrand. The l/r singularity of Z to be investigated in more detail in the next subsection is neutralize d by the measure dr” = 2r’dr’, ensuring an everywhere finite integrand. 8.
IMPROVEMENT OF THE GREEN’S FUNCTION
CONVERGENCE
FOR
THE
In spite of the fact that the above elimination of the direct beam now makes the integrals [Q. 1- .3] convergent even without the damping of the beam profile, the speed of convergence is nevertheless still unsatisfactory for numerical evaluation. In the integral expression [Q. 11, for example,_the factor of r~remaining from the measure d2v cancels the falloff of Q = Zr at T,J+~ and leaves only the 7 - ‘/’ falloff of the oscillatory Jo. The wasp-waist approximation, however, which eliminates the slow exponential damping, permits the following rewriting of the Hankel transforms, using Bessel’s equation and integration by parts, to increase the speed of convergence and thus make a rapid numerical evaluation feasible. 8.1. Zmproved Convergence of the Jo-Transform From the Bessel equation of order zero
d2Jo( rl) + 1 dJ,(rl) + Jo{11)= 0, h2
rl
(8-l)
ds
we may write
J&l) = -
d%( rl) 1 dJo(v) 7
-
-
11
Setting o( q/ r) = f(v) and inserting Equation 27rr2Q(r)
= - /j+)
d&(rl)-
-
drl -
(8.2) \ I
.
[8.2] into [Q.l]
somf(
yields
v)dJo( rl).
(8.3)
A partial integration converts the first term on the right into
-I;(q)qf(q)lo+~m;l(q){f(n)+?f.(q (8.4)
125
Forward Scattering Amplitudes
The first term vanishes as long as f(0) is finite and rlf(q) remains bounded as 7 + 03. The first term of the integral cancels the last integral of the equation above, leaving
27rr2Q(
4 = /j.f’(
7) 4o( 7)
= lo(4)vf+)lrn-
0
/“lo(v)
+Lf(?)).
0
P-5)
Assuming rlfl(q) is bounded as well, as will turn out to be the case in the forthcoming discussion, the first term again vanishes, leaving
t 0.1’1
8.2. Calculation of the b-radiance To calculate the scattered irradiance to insert into [Q.l’j is
I#-, zr), the appropriate function
(8.6) where C= F*exp( - u,zf),
w(q)=fi(7j/r,zf)=
VG arui
(8.7)
erf(s).
Now
f(v) = C*exp{w(s)}w’(rl)
(8.8)
126
JOHN LANGERHOLC
and
J,‘;fexp( - (~z’,~cx~)~) dz’
w’( 7) = @f
= uf zfexp ( - (w/242)
(8.9)
- o( rl)
so that
X
(-
(~/24~))
(8.10)
Hence the Hankel transform for the scattered irradiance is
X
-qw’(q)“+~‘(q)+$exp-(qn/2ar)~
dq.
(8.11)
8.3. Asymptotic Behavior of the k-radiance at Small r The most rapidly increasing contribution to the scattering n-radiance Zf(r, z) as r + 0 comes from the second term in the exponential expansion: ‘;( v,z)
= C( en(q~a)- 1) = Cj$
Q(v, z)m/m!,
with C= F*exp( - u,z).
(8.12)
Forward Scattering Amplitudes Inserting only the expansion,
127
m = 1 term
of this series into the inverse Bessel
=(oiC12~)SIdz’Smexp(-(llt’/2ar)‘)rllo(rg)d~. 0
(8.13)
0
and exchanging orders as has been done, we may evaluate the q-integral as a Laplace transform after a substitution of variables t =
?j2
p = ( 2’/2 Cq2 a=
r2/4
giving
00 =
J0 =
exp( - pt)Jo(2Jat)
(2ol/z’)‘exp(
dt =
p-le-‘lp
- (W/T?)“),
(8.14)
hence
I,(r, 2) = (Cu__cx2/7r)J,‘dz’exp(
-(rcy/z’)2)/z’2.
= (Cu__cz/7rr)[r,zdxexp(
= sexp(
- u,z)erfc(ar/z).
- x2)
(8.15)
JOHN LANGERHOLC
128
Thus the single-ray Green’s function has an asymptotic behavior proportional to 1 /T near the axis. 8.4. Calculation of the Tilt The first moment of the angular distribution
is calculated
by inserting
i)(S,Zf)=-if,(o,q,Z)i(‘1,Zf(Y=03),
(8.16)
where
f,(OA Zf) = - y-
(1-exp(
-(~.z@Y)~))>
(8.17)
into [Q.Z] with n = 1: M,,,[ I] (r, Zf) = bw(
?-, 4
=- J 1
27rr2
=&7(a/r)--;i?O
aJ( 7)
(8.18)
where
9(rl)=~(n)/(-iq,/q)=-~{1-rxp( -(~)“))f(vpzf) (8.19) 9(rllr)
Fuf r = - -exp(rl
2
a,zf)
e4rl)
.ii)
.
Inserting above yields
(Wf
with
(r, Zf.) =
(8.20)
129
Forward Scattering Amplitudes
Inserting Jl(s) dv = - 4drl) and integrating by parts improves the
convergence
of the expression:
x
d(n)(l--P{ - (Wf/2ar)2))
{
vf2 +-
2a2r2
exp
1( -
Wfl
2
ml-
)I) 2
(8.22)
8.5. Alternate Expression for ltilt If instead the limiting value 1 is subtracted from the exponential e“ to improve convergence, then there will be an extra term to Itilt, obtained by setting e* = 1, and equal to
J,-{1 -exp(
- (lz~/2ur)2]}lr(~)
ds
(8.23)
From standard tables of Hankel transforms, the last integral is found to have the value l-exp( where a=(+-/cxr)
‘tilt(
-l/a),
2. Hence the alternate expression for evaluating Itilt is
r,=exp(- (arlzf)2) +l-(I-exp(
-(‘lzf/2orr)‘))(eO(*)-I)J,(rl)dtl
(8.24)
An interesting feature of this approximation is that Z,,(O) = 1, hence the first E.-moment has an asymptotic behavior as l/r at small r, matching that of the scattering irradiance. As a consequence, this implies that (4,)
130
JOHN LANGERHOLC
will not approach 0 as r + 0, but rather lim (4,) = L r+O
(8.25)
cwJ?r.
Calculation of the Second Moments For the second moments, the quantities to be inserted into the inverse transform are given by 8.6.
Qz,o(l))=f2,o(ll)j(q,~fl~=0o),
(8.26)
Qo.,(~)=fo,~(‘l)j(tl,zfl~=~),
(8.27)
where
- (?(I-exp(
-(nr,i20)2))iL. (8.28)
fo,e(
11) =
-
f,(9
2P2
-(
q2-
7x2
1
=-+
tly
-“(O, 2 ffsqs
F(l-exp(
f,,
y(o,
7,
q)
8, x) + $$exp(
- (qz//2ar)‘]
-(qz~/2*)2]])2.
These each decompose into three integrals of different character. The constant terms l/2/3’ merely reproduce a multiple of the irradiance and need not be calculated separately. The remaining two are of the types [ Q.21 and [ Q.31 ,involving 119and T$, respectively. The formulas for the angular spread may then be written
(Ab)2=l/~z+{M,,(r)+ M&-)}/i(r-,q)
(8.29)
131
Forward Scattering Amplitds where
i%-,(v) =
2 z ( )I
ufzf
(y2
ew(- (vq/2a)2)
- y{l-exp(
-(qzf/2u)2])2
(8.30)
f(r,zf) I
(8.31)
(8.32)
2 OfZf
( ii
exp( - (wf/2~)2}
-f-$(1-exp{
-(?JZf/201)2]]2
A&,(rl)= :
(112
i(r,zf).
(8.33)
I
The integrals with
involving
7: are most conveniently evaluated using [Q.3]
(8.34)
rlzf ,2a)2])2]
i(r,
zf),
(8.35)
132
JOHN LANGERHOLC Fexp(
9ee( V/T) = c
-(tlzf/2cYq2)
i -~{~-~~~{-(n~f,2~~)2)}2
eucs),
(8.36)
I (8.37)
we find
A4&-)=(C/2*r%2 'J - 2u$zr2
0
O1 afZfexp (-(l)zf/2Ur)2]
(8.38)
i
[I-exp(
-(qzf/2ar)2]]2
e”(@J,(v)dv. I
In the second integral the factors of J1 fall off either exponentially fast as tm2. To evaluate MO,, we apply [Q-2] with n = 2 and
9&)
= - ‘(z2zf)
i( r, zf)
or as
(8.39)
and use the fact that dJ,(q)/dq = - J1(q) to obtain
M&r)
= (C/2rr2az
)[av~(n)eu(v)41(n).
(8.40)
Sufficient falloff can be obtained by integrating once by parts
Mar(r)
= - (C/2rr2a2
) L=Jb) dq
~{~(1))+r)~‘(rl)+rl~‘(rl)~(rl)}e~(”). (8.41)
Forward Scattering Amplitudes
133
To use [Q.2] to evaluate M,,,
4rr(9) = -
for which n = 2 and
exp{ - (rIz//2o)‘)
y
1 - 7
{ I-exp(
- (qzf/20)z]]2)
i(r,
zf),
(8.42)
we rewrite, from Bessel’s equation (8.43)
-~on(rl)=Jo(rl)+J~(rl)lrl=lo(9)-J~(rl)lrl. Insertion yields, finally,
M,,(r)
= ( C/27rr2a2
x
urzfexp [ - (11+q2] i
- 2u-J2r2
(1-exp{
eWcrl). (8.44)
-(~zf/2ar)2])2
1 9.
RESULTS
The Gaussian parameters were computed as functions of the radial distance r from the beam axis in the focal plane zf = 2 km for a laser of F = 1 megawatt for several values of scattering parameters in the wasp-waist approximation. The maximal irradiance values were estimated by the Gaussian fit and plotted in Figure 3 for the values ur = 0.035 km-’ and o = 15 radd’. Owing to lack of data on the absorption coefficient, the extinction was arbitrarily set at ue = 2 a, in all computations. The values can be adjusted for different values of the absorption coefficient by multiplying with the appropriate exponential. It is seen that the maximal radiance value is exceeded within 4 meters of the beam axis. Figure 4 shows a plot of the single-ray Green’s function for the irradiante for the values uf = 0.035 km-’ and (Y= 8.11 rad-‘. Figure 5 shows the same function for the larger value of uf = O.l6/km. In spite of the
JOHN LANGERHOLC
134
(W/m%) 4
lOe?--
lOe6-R d" ' a n c e
lOe5--
10a4--
10e3
! 0
b
I
2
4
6
6
10
12
14
16
16
20
Radial Distancefmm Axis
FIG. 3.
(Wlm2)
I
0
50
l (m) 100
150
200
250
Radial Distance from Axis
FIG. 4.
300
350
Cm)
135
Forum-d Scattering Amplitudes (Wlmz) 1000.00
1 I
r
r
100,oo --
d”
I
a n
IO,00 --
C e
I,00 3
0
5
10
15
20
25
30
35
40
45
50
ä(m)
Radial Distance from Axis FIG. 5.
increased absorption coefficient the values of Z are approximately three times as great. The angular spreads A+ of the distributions are generally small. They are less than 0.1 rad at r = 10 m for both values of the scattering coefficient u and in both r and 8 directions. performed numerically in the The power integr ai s 2 rr / Z(r, zf)rdr focal plane up to r = 370 m agree with the theoretical values integrated over the entire focal plane to within 4% in both cases. The author is indebted to G. Born and H. Hermansdorfer of MesserSchmitt-Bijlkow-Blohm, 8012 Ottobrunn, West Germany, for financial support of this work, and especially to G. Born for suggesting the problem and its relevance to questions of laser safety. He wishes to thank M. Gorriz, also of MBB, for discussions and assistance with the computations, and P. Glijckl for the figures.
APPENDIX.
LIST
OF
SYMBOLS
Z(+,r, z): The radiance
as function
of the angular orientation
coordinates (r, 2). Z(r, z): The irradiance, the spatial coordinates
or angular (r, z).
integral
of the radiance,
4 and spatial as function
of
JOHN LANGERHOLC
136 f(f, 1, z):
The
Fourier
transform
of the
radiance
in angular
and radial
variables.
f(q, z): The Fourier = i(O,q, 2). CY:Inverse breadth
transform
of the irradiance
of the forward
in the radial
variables
peak of the phase function.
ur: The coefficient of forward scattering. (Y,: Inverse breadth of the reverse peak of the phase function. a,.: The coefficient of backscattering. a,,: The coefficient of isotropic Us: The absorption coefficient. a,: The extinction coefficient. p: Inverse y: Inverse
scattering.
angular divergence of the input beam. beam radius RfoL in the focal plane z = zf.
Rfoc: Beam radius (l/e) in the focal plane x = zf. Ri,: Beam radius (l/e) in the input plane. zf: Distance from the input plane to the focal plane.
REFERENCES 1 2 3
4 5 6 7 8
9
10 11
H. S. Snyder and W. T. Scott, Phys. Reo. 76:220-225 (1949). L. R. Dolin, Propagation of a narrow beam of light in a medium with strong anisotropic scattering, Zzv. Vyssh. Uchebn. Zaved. Radiofiz. 9:61- (1966). R. L. Fante, Propagation of electromagnetic waves through turbulent plasma using transport theory, ZEEE Trans. Antennas Propagat. AP-21:750-755 (1973). A. Ishimaru, Wave Propagation and Scattering in Random Media, Academic Press, New York, San Francisco, London, 1978. A. Zardecki and W. G. Tam, Laser beam propagation in particulate media, Can. 1. Physics 57:1301- (1979). W. G. Tam and A. Zardecki, /. Opt. Sot. Am. 69:68-70 (1979). W. G. Tam, Multiple scattering corrections for atmospheric aerosol extinction measurements, Appl. Optics 19:2090-2092 (1980). W. G. Tam and A. Zardecki, Multiple scattering corrections to the Beer-Lambert law. 1: open detector Applied Optics 21:2405-2412 (1982); 2: detector with a variable field of view, A&. Optics 21:2413-2420 (1982). A. Deepak, V. 0. Farrukh, and A. Zardecki, Significance of higher-order multiple scattering for laser beam propagation through hazes, fogs, and clouds, Applied Optics 21:439-447 (1982). J. Langerholc, A beam-width formula derived from a two-component scattering theory, in preparation. J. Langerholc, Volumes of diced hyperspheres: resumming the Tam-Zardecki formula, Appl. Math. Comput. 30:1-18 (1989).
Forward 12
13
14
Scattering
Amplitudes
137
‘K. Altmann, Forward-scattering formula of Tam and Zardecki evaluated by use of cubic sections of spherical hypersurfaces, Appl. Optics 28:4077-4087 (1989). J. A. Weinman, J. T. Twitty, S. R. Browning, and B. M. Herman, Derivation of phase functions from multiply scattered sunlight transmitted through a hazy atmosphere, J. Atmos. Sci. 32:577- (1975). A. Zardecki and S. A. Gerstl, Multi-Gaussian phase function model for off-axis laser beam scattering, Appl. Optics 26:3000-3004 (1987).
‘Although perhaps now only of academic interest, this paper could have been considerably shortened and the unnecessary approximations avoided by using the original recursion relation of [ll] to compute the volume functions iteratively. In fact, the area functions themselves need never be computed if a Stieltjes integral is used, since the appropriate measure Surf.(r) dr = dVol.( r) is in fact the differential of the volume function.