Arizona Center for Mathematical Sciences, University of Arizona, Tucson, Arizona 85721
Also with the Optical Sciences Center, University of Arizona,
Tucson, AZ 85721
Also with University of Ostrava, Ostrava, Czech Republic
M. Mlejnek 
, E.M. Wright , and J.V. Moloney
We present numerical simulations of nonlinear pulse propagation in air whereby an initial pulse is formed, absorbed by plasma generation, and subsequently replenished by power from the trailing edge of the pulse. This process can occur more than once for high power input pulses and produce the illusion of long distance propagation of one self-guided pulse.
The observation of long distance apparently self-guided pulses in air [1, 2, 3, 4, 5] has attracted much recent interest, with applications in communications and lightning channeling envisioned [3]. The basic phenomenon characterized by long filaments and propagation over tens of meters is not in dispute but the physical mechanism underlying it has still to be uncovered in detail. Early discussions centered on the idea that the nonlinear self-focusing (SF) in air would be countered by the defocusing effect of the electron plasma generated by multi-photon absorption (MPA), and result in a stable self-guided beam. Numerical simulations showing stabilization of the peak field intensity were presented in Refs. [1, 6], and Brodeur et. al. [5] have used the moving focus model to explain the features of their experiments.
In this Letter we present numerical simulations of nonlinear pulse propagation in air to elucidate the physical mechanisms involved. Our study includes the effects of SF, MPA, group velocity dispersion (GVD), and absorption and defocusing due to the generated electron density [7, 8]. Although our simulations have not been validated against current experiments, it was clear that they do not yield results resembling the self-guiding picture. Instead a dynamic picture emerges in which pulses form, are absorbed, and are subsequently replenished by new pulses, thereby creating the illusion of one pulse which is self-guided.
Our model for pulse propagation in air is a nonlinear extension of
one due to Feit and Fleck [7].
Assuming propagation along
the z-axis, the equation for the slowly varying electric field envelope
in a reference frame moving at the group velocity is
[8]:
where the terms on the right-hand-side describe transverse
diffraction, GVD, absorption
and defocusing due to the electron plasma, MPA, and nonlinear SF.
Here 
is the optical frequency, 
the intensity, 
,
,
is the electron density,
the cross-section for inverse
bremsstrahlung, 
is the electron collision time,
is the K-photon absorption coefficient, and
the nonlinear
change in refractive-index for a continuous wave (cw) field
is 
. The critical power for self-focusing
collapse for cw fields is then 
.
The normalized response function
R(t) accounts for delayed nonlinear effects, and f is the
fraction of the cw nonlinear optical response which has its origin in
the delayed component.
The evolution of the electron density is described by the Drude model
[7, 8]
The first term on the right-hand-side
of this equation describes growth of the electron plasma by cascade
(avalanche) ionization, the second term is the contribution of MPA,
and the third term describes radiative electron recombination.
We are interested in solutions of these equations for an input
collimated Gaussian beam
where 
is the peak input power,
the spot size,
characterizes the pulse length,
and 
is the Rayleigh range of the input beam.
For the numerical simulations presented here we chose an operating
wavelength 
nm, 
mm,
for which the Rayleigh range is 
m, and 
fs.
The material parameters employed in our simulations
for air at standard temperature and pressure are
as follows: 
, 
m 
/W
[9], giving a critical power of 
GW,
eV representative of the constituents of the air,
which yields K=7 for the order of the MPA, and
m 
W 
for the MPA
coefficient from the Keldysh theory of MPA [10],
s (typical value from [7]),
m 
/s,
m was calculated using
[7].
The GVD for air was taken as
k''=2.0 fs cm 
and we shall comment on using other values.
The numerical results were checked by doubling the space
and time resolutions which led to no significant changes in the
behaviour of the results presented here.
Nibbering et. al.
[9]
have shown that air
exhibits a delayed nonlinear response due to rotational Raman
scattering with f=1/2. Their work shows that
the response function R(t) initially undergoes a damped harmonic
motion, followed by spontaneous resurgences for times beyond 2 ps.
For the sub-picosecond pulses considered here the resurgences are
not relevant and we model the response function using a damped
oscillator model
,
where 
, and we have used
THz
and 
THz to mimic the exact response function given in Refs.
[9].
Fig. 1 shows the global maximum over time
of the on-axis intensity as a function of propagation distance z for
peak input powers of 5.5 
(dash-dashed curve),
6.0 (long-dashed curve),
and 6.5 (solid curve).
In each case the maximum intensity initially grows explosively due to SF,
but is then limited by MPA and absorption and defocusing due to the
electron density [1, 6].
After limiting, the intensity remains fairly constant up
to a distance of around 2 m, or roughly a Rayleigh range, after which it
decays.
This observation is in agreement with the moving focus prediction by
Brodeur et. al. [5] that light filaments should end at a
Rayleigh range of the input beam, an exception being the case
which we will comment on later.
Fig. 1 reveals a notch in the curves where
dI/dz is discontinuous.
To analyze this we show the on-axis pulse profiles in Fig. 2(a) for
z=0 m (dotted curve), z=0.55 m (dashed curve), and z=1.1 m
(solid curve), for the case , and
we see that the temporal compression effect due to SF
[11]
has contracted the pulse resulting in a single peaked pulse.
At this stage the peak intensity growth is
limited but the temporal pulse compression still proceeds until it is stopped
by the normal GVD
[12, 13].
In general, if the normal GVD is too small this compression
proceeds until the numerical scheme breaks down and/or the envelope
approximations underlying the model break down.
We plan to investigate other linear and/or nonlinear terms which
may become important in this regime, e. g. shock and nonparaxial
terms, in the future.
For the value of GVD employed here the pulse compression is arrested
before the model and numerics are violated.
For a propagation distance of z=1.1 m (solid curve) we see that the
pulse profile shows two peaks separated temporally by
fs.
The leading peak occuring at earlier times is directly associated with
the single peak at z=0.55 m, and the trailing peak developed out of
the background.
Upon further propagation the leading peak decays while the trailing
one remains, until it too decays at around a Rayleigh range.
The notch in the intensity curve is therefore explained as that
propagation distance where the increasing trailing peak takes over as
the global maximum from the decreasing leading peak.
We also note that a similar notch was experimentally and
numerically observed in the filament energy measurements of ref.
[5].
The double pulse profile in Fig. 2(a) is reminiscent of the pulse-splitting
phenomenon for nonlinear pulse propagation
with normal GVD
[12, 13].
Here normal GVD plays a role in that it arrests the temporal
compression,
but the leading pulse develops first
followed by the trailing pulse, whereas in the pulse-splitting case the
two pulses develop simultaneously.
In analogy to pulse-splitting the
signature of the double pulse in Fig. 2(a) will be a modulated spectrum
[13],
with period 
, or for T=200 fs,
nm. Fig. 2(b) shows the spectrum of the double pulse
in Fig. 2(a)
and we see the development of a spectral modulation with 10 nm period.
To expose the physics of the development of the trailing edge pulse,
we note that the leading edge pulse at z=0.55 m
shown in Fig. 2(a) (dashed curve) has its maximum displaced by -50 fs
from the input pulse (dotted curve).
This occurs because as the input pulse
focuses and the intensity increases, an electron density is generated via
MPA and cascade ionization, which produces a defocusing effect which is
more pronounced for the trailing portion of the pulse.
This is confirmed by the fact that the trailing edge of the pulse
(dashed curve) is reduced with respect to the input (dotted curve).
To elucidate further Fig. 3 shows the spatial profiles for propagation
distances 
m and for the time slice t=70 fs, which is
that at which the trailing pulse forms (see Fig. 2(a)). The input
beam profile is Gaussian (dotted curve), but for z=0.55 m (dashed curve)
the central intensity is suppressed with respect to the input, the
beam power having been displaced into a spatial ring by the defocusing
effect.
However, for z=1.1 m (solid curve) the beam
profile has relocalized on-axis:
This spatial defocusing of the trailing portion of the pulse into a
ring and subsequent relocalization into a strong peak via SF is the
physical mechanism by which the trailing edge pulse forms.
Refocusing of the trailing edge of the pulse was also noted by
Kosareva et al. [14] (they do not include GVD) using similar
model including tunneling ionization.
There are also indications from
experiments that the energy outside the self-guided channel
is important for the prolonged propagation of pulses in air
[4].
The process of pulse decay and replenishment can occur more than once:
The solid curve in
Fig. 1(a) for
shows the growth of a second trailing
edge pulse beyond the Rayleigh range m, and this is again due
to the relocalization of power that is initially displaced into a
spatial ring by the defocusing effect of the plasma.
We remark that the occurence
of the second trailing edge pulse in our numerical simulations
was dependent upon including the delayed nonlinear response due to
the rotational Raman scattering.
In summary, we have presented numerical simulations to expose the physics underlying long distance propagation in air. We have shown that although the maximum intensity stabilizes over roughly a Rayleigh range, in agreement with the moving focus model of Brodeur et. al. [5], the evolution of the pulse is very dynamic, involving the development of a leading edge pulse which subsequently decays and is replaced by a new pulse. In addition, our simulations indicate that the inclusion of rotational Raman scattering [9] permits high power pulses to continue the process of pulse decay and replenishment beyond the Rayleigh range limit. There is also a potential link between the spatial ring formation predicted in the present work and the conical emission observed in recent experiments [4].
Acknowledgement: Effort sponsored by the Air Force Office of Scientific Research, Air Force Materiel Command, USAF, under grant numbers AFOSR-97-1-0002 and AFOSR-94-1-0463. Ewan Wright was also supported by AFOSR contract F49620-94-1-0343.
