1. Introduction

The current development trend for short pulse laser systems is towards higher repetition rates (>kHz) and average power (> kW). The high average power, high peak power pulse trains produced by these systems drive a unique physics regime of atmospheric propagation and material interactions that has outpaced predictive modeling capabilities. Specifically, the ultrashort pulse trains possess average powers near those of high-energy lasers (HELs) and peak powers sufficient for nonlinear propagation. This high-average and high-peak power regime results in both linear and nonlinear heating of the air, combining, for the first time, high-average-power effects like thermal blooming with those found in filamentation, e.g. ionization and rotational excitation [1–6]. The propagation of laser filaments in this new regime has not been thoroughly investigated.

Previous efforts have focused either on low peak power HEL propagation where traditional thermal blooming plays a role [7], single pulse, high peak power pulse propagation such as filamentation studies [3,8–11], or experiments with high peak power trains with relatively low average power ($\approx$ 1 W) [5,6]. Here, we present a theoretical framework suitable for implementation in a propagation code to model high average power ( kW), high peak power pulse trains.

We consider a train of pulses, each of spot size $R_0$, duration $\tau _p$, and peak power $P_0$, separated by duration $\tau _s$. Each pulse deposits thermal energy into the air which leads to a small evacuation of the air density on axis and a consequent decrease in the refractive index. Between pulses, the air cools and the refractive index increases towards its equilibrium value until the next pulse comes through. If $\tau _s$ is comparable to or shorter than the cooling time (either through convection or conduction), the refractive index perturbation can build up over the pulse train. This is shown schematically in Fig. 1 along with representative values of the pulse duration, $\tau _p$, peak power, $P_0$, separation time, $\tau _s$, acoustic transit time, $\tau _a$, and thermal conduction time, $\tau _D$.

 

Fig. 1. Schematic diagram of a laser pulse train propagating to the left. Solid curves represent laser power, the dashed curve represents the refractive index of air.

Download Full Size | PDF

2. Absorption of laser energy by air for ultrashort laser pulse

Laser pulses propagating in air lose their energy to linear absorption and scattering, ionization, rotational and vibrational Raman excitation of air molecules, and for pulses longer than the electron-neutral or electron-ion collision time, inverse Bremsstrahlung. All of these mechanisms, except scattering, result in heat being deposited along the propagation path of the laser. When an air molecule is ionized a plasma is created and the laser pulse loses energy equal to the ionization energy, $U_I$, of the molecule plus, in the case of above-threshold ionization (ATI), whatever residual energy the free electron carries away. If the electron remains in the laser field and collides with a neutral or ion, it gains approximately the oscillation energy. After the laser pulse has passed, energized electrons transfer their energy to air molecules on the collisional time-scale of several tens of picoseconds. On a longer nanosecond time-scale, electrons recombine with and attach to ions and neutral molecules, respectively. For a weakly-ionized plasma at atmospheric pressure, the primary electron loss mechanisms in air are non-radiative, exothermic, three-body processes [12] which lead to local air heating. The laser pulse can also directly excite rotational and vibrational modes whose energy eventually becomes heat on a molecular collision time-scale.

The depletion of energy from the laser pulse can be described by the equation $\partial F_L / \partial z = -\alpha F_L$, where $F_L$ is the laser pulse fluence and the effective extinction coefficient, neglecting scattering losses, can be written as $\alpha = \alpha _0 + \alpha _{ion} + \alpha _p + \alpha _R$ where, in order, the terms denote the contributions from linear absorption, ionization, plasma heating (inverse Bremsstrahlung), and rotational Raman excitation. The change in laser fluence with $z$ is also equal to the negative of the energy density deposited by the laser into the air, which we write as $u_a = u_0 + u_{ion} + u_{p} + u_R$ where the various terms on the RHS are the energy density deposited into the various process denoted by the subscripts, as in the expression for $\alpha$. With this, the absorption coefficient for the process labeled $j$ is given by $\alpha _j = -u_j/F_L$. For nonlinear absorption processes, $\alpha$ is a function of the local laser intensity. One can also define a macroscopic absorption coefficient $\Lambda$ by integrating the equation for fluence loss over transverse coordinates to give $\partial U_L / \partial z = -\Lambda U_L$, where $U_L = \int F_L d^2r$ is the energy of the laser pulse, and $\Lambda = U_L^{-1}\int \alpha F_L d^2r$. In the following sub-sections, we derive simple expressions for the various contributions to the laser absorption coefficient that can be used in propagation codes that model the interaction of many pulses in a pulse train.

2.1 Ionization and electron heating

As the laser pulse creates and energizes free electrons, it loses fluence according to

We now evaluate $n_e$ and $u_e$. The equation governing the density of free electrons, for pulses short compared to the electron recombination and attachment times is $\partial n_e/\partial t = \gamma _{ion}(I)(n_n-n_e)$, where $n_n$ is the ambient neutral density, $\gamma _{ion}$ is the photoionization rate, and we have neglected collisional ionization on the short time-scales of the laser pulse. The evolution of the electron energy density, $u_e = 3 n_e k_B T_e/2$, is given by $\partial u_e/\partial t = \left \langle\vec {J}\cdot \vec {E}\right \rangle$, where $\vec {E}$ is the laser electric field, $\vec {J}$ is the linear part of the electron current governed by $\partial \vec {J}/\partial t + \nu _e \vec {J} = (\omega _{pe}^2/4\pi )\vec {E}$, $T_e$ is the electron temperature, $\omega _{pe}^2 = 4\pi e^2 n_e/m_e$ is the plasma frequency squared, $\nu _e$ is the electron collision frequency, and $k_B$ is Boltzmann’s constant. Here, we have assumed that we can neglect electron recombination and cooling by interaction with the air over the duration of the laser pulse.

We can write a simple, closed form solution for the electron energization by assuming that $\nu _e$ and the laser intensity are constant over the duration of the pulse, and that the pulse has a simple radial profile, i.e., $I(r,\tau ) = I_0 f(r)$ for $0<\tau <\tau _p$ and zero for all other $\tau$. In this case, the electron number density and energy density immediately following a laser pulse of duration $\tau _p$ are

The absorption coefficients due to ionization and electron heating are given by the respective right-hand terms of Eq. (1) divided by $-F_L$, i.e., $\alpha _{ion} = n_0 U_I n_e/F_L$ and $\alpha _p = u_e/F_L$, with $n_e$ and $u_e$ given by Eqs. (2) and (3). Assuming a Gaussian radial profile, $f(r) = \exp (-2 r^2/R_0^2)$, the macroscopic absorption coefficients due to ionization and plasma heating are given by

2.2 Rotational Raman excitation

Laser energy loss due to rotational Raman excitation follows from the slow variation in time of the nonlinear Raman polarization field. It can be derived from the paraxial wave equation for the laser envelope $\vec {A}$ [4,16], shown here in a simplified form only containing self-phase modulation effects

The equation for laser fluence loss to rotational Raman excitation derivable from Eq. (7) is given by

2.3 Comparison of absorption processes

We compare the various magnitudes of the absorption coefficients derived previously by plotting the absorption coefficients as a function of pulse length for a fixed pulse energy. Figure 2 shows the results for a 3 mJ pulse, and a 10 mJ pulse with a focused spot size of $R_0 = 200 \mu$m. For the 3 mJ pulse, rotational Raman is the dominant absorption mechanism for pulse durations $>200$ fsec, while ionization loss is dominant for pulse durations $<100$ fsec. Absorption due to electron heating is negligible over the range plotted. For the 10 mJ pulse, ionization loss dominates for pulse durations $<0.3$ psec and Raman excitation for pulse durations $>1$ psec. Absorption from electron heating is only plotted for pulse durations greater than the electron collision time, and is seen to be comparable with ionization loss. Linear absorption due to air molecules is indicated by the dashed black line. The lower panels plot the electron temperature versus pulse duration and show that the electrons are energized to approximately 1 eV for the 3 mJ pulse and several eV for the 10 mJ pulse. Making use of the expression $u_e = 3n_ek_BT_e/2$ and Eqs. (2) and (3) it follows that

 

Fig. 2. Electron density (solid black curves), and macroscopic absorption coefficients, $\Lambda _i$, due to Raman excitation (blue curves), ionization (green curves), electron heating (red curves), and air molecules (black, dashed line, $\alpha _0=1.6\times 10^{-4}$ km$^{-1}$ at $\lambda _0=800$ nm) versus pulse duration for a fixed pulse energy, a) 3 mJ, and c) 10 mJ. Electron temperature vs. pulse duration is also plotted for pulse energies 3 mJ (panel c), and 10 mJ (panel d). For all figures $R_0 = 0.2$ mm.

Download Full Size | PDF

3. Air heating and thermal refractive index perturbation caused by an ultrashort pulse

Some of the energy lost by the laser pulse contributes to heating air molecules in the laser path, resulting in a refractive index perturbation that can affect subsequent pulses. Physically, the laser pulse produces a region of elevated temperature that is initially without a density perturbation. An acoustic wave is then generated that carries the excess energy away from the heated region and leaves behind a heated region in pressure balance with a density depletion which evolves on a hydrodynamic time scale [6,10,18]. In this Section we calculate the refractive index perturbation left behind by one pulse over a time long compared with the acoustic transit time across the propagation path, but comparable to the convection and thermal conduction times.

The change in refractive index of a gas due to a density perturbation $\delta \rho$ is given by $\delta n = (n_0 -1)(\delta \rho /\rho _0)$, where $n_0$ is the ambient refractive index of the propagation medium, and $\rho _0$ is the ambient mass density. In general, when the density and temperature perturbations caused by the laser in a gaseous medium are small, they are governed by the linearized fluid equations

There are two parameter regimes that can be considered, which are dictated by timescales that are greater or less than $\tau _a = R_0/c_s$, i.e., the acoustic transit time across the pulse. For time scales short compared with $\tau _a$, the propagation of the acoustic wave is important, convection and conduction can be neglected, and Eqs. (16) and (17) reduce to

For the purpose of analysis, we assume, as before, that a pulse has a flattop longitudinal intensity profile, $I(r, \tau ) = I_0 f(r)$ for $\tau < \tau _p$ and zero otherwise, where $\tau =0$ marks the arrival of the pulse, and $f(r) = \exp (-2r^2/R_0^2)$. Since $\tau _p << \tau _a$, we can neglect the operator $c_s^2\nabla ^2$ in Eq. (18) for times comparable to $\tau _p$ and integrate in time to yield, $\partial _t \delta \rho = (\gamma -1)\alpha I_0 \nabla ^2 f(r)t^2/2$. Substituting this result into Eq. (19) and integrating in time gives

For $\tau _p < t << \tau _D$, we can similarly solve $\left (\partial _t^2 - c_s^2 \nabla ^2\right )\partial _t \delta \rho = 0$ for $\partial _t \delta \rho$, using $\partial _t \delta \rho = (\gamma -1)\alpha I_0 \nabla ^2 f(r)\tau _p^2/2$ as the initial condition, then substitute into $\partial _t \left (\delta T/T_0\right ) = \left (\gamma - 1\right )\partial _t \left (\delta \rho /\rho _0\right )$, and integrate to yield

4. Inter-pulse thermal interaction of a pulse train

Given the refractive-index imprint of a single pulse, we now consider a train of pulses with wavelength $\lambda = 800$ nm and duration $\tau _p < 1$ psec, with spot size $R_0 \approx 0.2$ mm and separation $\tau _s \approx 1$ msec ( kHz rep-rate) propagating through air at STP, which is characterized by $\alpha _0 = 1.6\times 10^{-4}$ km$^{-1}$, $\rho _0 = 1.23$ kg/m$^3$, $C_p = 1$ kJ/(kg K), $\gamma = 1.4$, $\chi = 21.7$ mm$^2$/s, $n_0 - 1= 2.9\times 10^{-4}$, and $T_0 = 293$ K. The acoustic speed is $c_s = 343$ m/s and $\tau _a \approx 0.6 \mu \text {s}$, which is much longer than $\tau _p$ so that heating of the air occurs in the $t^3$ regime. Between pulses, cooling will occur in the isobaric regime for times $\tau _a < t < \tau _s$. The characteristic time for thermal conduction of a Gaussian perturbation with a $1/e^2$ radius of $R_0$ is $\tau _D = R_0^2/(\chi ) \approx 2$ msec, which is comparable with the inter-pulse separation.

4.1 Convective vs diffusive cooling

Before proceeding, it is useful to estimate the relative effects of convective and conductive (diffusive) cooling of the laser-heated air between pulses. Neglecting transverse wind, the relative magnitudes of thermal convection (by buoyancy) and conduction are parameterized by the Peclet number, $P_e = v_B R_0/\chi$, where $v_B$ is the buoyancy velocity, and $P_e>>1$ describes the convection-dominated regime. The buoyancy velocity can be estimated from $(dv_b/dt) \approx (\delta \rho /\rho _0)g$, where $g$ is the gravitational acceleration, and the approximation $\delta \rho /\rho _0 << 1$ was made. This expression for the buoyancy velocity assumes that viscous stresses are negligible. The Reynolds number $R_e=2R_0v_B/\nu$, computed with this velocity and taking $\nu =0.15$ cm$^2/$s for the kinematic viscosity of air, is found to have a value of 3700, sufficiently large to justify the assumption.

If $\tau _s$ is much shorter than the convection or conduction times, then $\delta \rho /\rho _0 \approx (\delta \rho _1/\rho _0)(t/\tau _s)$, where $\delta \rho _1$ is the density change due to a single pulse. Letting $t = R_0/v_B$, i.e., the convection time, we obtain $v_B \approx [|\delta \rho _1| g R_0^2/(\rho _0 \tau _s)]^{1/3}$. From Eq. (23), we have $\delta \rho _1 \approx \alpha I_0 \tau _p /\left (C_V T_0\right )$ on-axis, after a few acoustic transit times. Using the results for $v_B$ and $\delta \rho _1$, the Peclet number is

Figure 3 plots $P_e$ vs pulse duration for a pulse train with fixed energy per pulse $U_L = 5$ mJ for various average power, obtained by varying the pulse separation time. For shorter pulses with increased absorption, convection due to buoyancy can be comparable to thermal conduction. Hence, for these shorter pulses, the heating of the air column should result in a deflection of the pulses downward as heat is convected upward. This will be seen in the simulations presented in the following sections.

 

Fig. 3. Peclet number versus pulse duration for a pulse train with fixed energy per pulse $U_L = 5$ mJ with average power 500 W (solid curve), 50 W (longer-dashed curve), and 5 W (shorter-dashed curve) for $R_0 = 200 \mu$m.

Download Full Size | PDF

4.2 Heating by a pulse train

For the purpose of obtaining a simple analysis, we consider the regime where $P_e < 1$ and the cooling between pulses is primarily due to thermal conduction. For pulse durations that are much shorter than the thermal conduction time, heating and cooling can be considered separately in time. If a single pulse produces a thermal index change $\delta n_1 \exp (-2 r^2/R_0^2)$ then the pulse immediately after it experiences an index $\delta n_1 (1+ \tau _s/\tau _D)^{-1}\exp \left [-2 r^2R_0^{-2}(1 + \tau _s/\tau _D)^{-1}\right ]$ [19], where we have assumed that the heated region diffuses transversely into free space.

Consider the change in refractive index on-axis. As an approximation, we assume that the transverse shape of the heated region under the influence of many pulses, for the purpose of estimating the effect of conduction, remains Gaussian near the axis with a spot size $R_0/\sqrt {2}$. Hence, the on-axis change in refractive index caused by the $N^\text {th}$ pulse and experienced by the pulse immediately following it (i.e. at a time $\tau _s$ after the $N^\text {th}$) is $\delta n_N = \delta n_{N-1}\left (1 + \tau _s/\tau _D\right )^{-1} + \delta n_1$. This recursion relation can be written

 

Fig. 4. Normalized refractive index change from Eq. (26) as a function of $N$ for $\tau _s/\tau _D = 0.1$ (solid curve), 1 (longer-dashed curve), and 5 (shorter dashed curve).

Download Full Size | PDF

5. Simulations

To simulate the propagation of a train of pulses undergoing filamentation, we used NRL’s PyCAP (Python Code for Atmospheric Propagation). The simulation solves each pulse in 3D using the Nonlinear Envelope Equation (NEE) of Brabec and Krausz [20], modified to include refractive index perturbations due to air heating, as well as plasma and ionization effects. The evolution of the electric field is given (in Gaussian units) by:

PyCAP solves the NEE using a symmetrized Fourier split-step with an adaptive propagation step which bounds local error [22]. For these simulations, the nonlinear polarization density envelope is $p_{NL} = p_{elec} + p_{rot} + p_{free} + p_{ion}$. The electronic Kerr term is given by $p_{elec}=(n_0/2\pi )n_2IA$, while the rotational Raman response $p_{rot}$ is described in detail in Eqs. (7–9). The polarization density due to free electrons is given by [23]

All energy lost by the pulse is assumed to heat the air, on a time scale short enough so that Eq. (23) can be used to approximate the induced refractive index change for each pulse, with the quantity $\alpha I_0 \tau _p$ in Eq. (23) replaced by the more general expression for lost fluence,

In between pulses, the perturbed index and velocity evolve according to the 2D linearized fluid equations in the isobaric limit,

Filamentation of a single pulse

Before attempting to simulate a train of pulses undergoing filamentation, we first verify that we can model the filamentation of a single femtosecond pulse. For comparison with experiment, we choose parameters identical to those found in Fig. 2(a) of [25], in which peak electron density was measured experimentally as a function of propagation distance. Pulses at $\lambda _0=800$ nm were focused at $f/240$ through a $f=95$ cm lens. The power was kept at a fixed 17 GW, while the pulse duration was varied from $\tau _p=40$ fs to $\tau _p=120$ fs. The simulation results are shown in Fig. 5 alongside the experimental measurements for the same parameters.

 

Fig. 5. (a) Experimental and (b) simulated peak electron density in a filament as a function of propagation distance, in comparison with Fig. 2(a) of [25].

Download Full Size | PDF

From Fig. 5, the simulated filament length and onset distance match the experiment to within about 5%. The maximum peak electron density matches very well in the 120 fs case and within 40% in the 40 fs case. In the 120 fs case, the secondary peak is absent from the simulation, which may be due to differences between simulation and experiment in the temporal pulse profiles, which have been shown to significantly affect ionization and energy deposition [26].

Simulation of kHz filamentation

To demonstrate the validity of our models and as an example of how the thermal interaction of a pulse train can lead to observable propagation effects, we again model parameters from a previously-published experimental work [5]. The authors focused pulses to filament in air (at 1 atm) and Xenon (at 1 atm and in a gas cell at 2.7 atm) at 1 kHz. They found that the heated air channel thus created rises due to buoyancy as it simultaneously expands due to thermal diffusion. This eventually creates a refractive index perturbation which pulses experience as a diverging lens offset transversely from the propagation axis. As the index perturbation builds up from the repeated filamentation, pulses are increasingly deflected downward as they pass through the heated region.

As in [5], we use 45 fs FWHM, 1.5 mJ pulses, focused at $f=60$ cm through a 4mm diameter aperture, with temporal and spatial profiles as described previously, and operating at a repetition rate of 1 kHz. Several quantities of interest are shown in Fig. 6. In particular, the electron density profile has a longitudinal FWHM of approximately 9 cm, which we take to be the interaction length for the purpose of comparison with interferometric measurements.

 

Fig. 6. (a) Peak electron density, (b) plasma channel FWHM, (c) laser fluence FWHM, (d) energy deposited, and (e) the peak intensity, as a function of propagation distance for the first (solid line), second (dash-dotted), third (dotted) and fourth (dashed) pulses, using the parameters from [5].

Download Full Size | PDF

The first pulse is notably different from subsequent pulses in all of the quantities shown in Fig. 6, in that the filament is significantly shorter longitudinally than for the pulses which follow. This is presumably due to the approximately axially-symmetric density perturbation as described in Section 4, which acts as a defocusing lens, in opposition to the effects of Kerr self-focusing and the phase induced by the 60cm lens at the aperture. This axisymmetric thermal lens reaches a steady-state by approximately the fourth pulse. From Fig. 6, in the area of peak heating between 7cm and 4 cm before the geometric focus, the fluence FWHM is $\approx 150 \mu$m, corresponding to $\tau _s/\tau _D\approx 4$. Hence, the number of pulses required for the thermal lens to reach a steady-state is in good agreement with the theoretical results shown in Fig. 4.

After examining the thermal interaction of the first few pulses, we then propagated 500 pulses through filamentation, with each pulse depositing energy on the 2D refractive index screens while simultaneously having its phase modified. Figure 7 shows the location of the transverse (fluence) centroid in the direction parallel to gravity (negative is down) for every fortieth pulse as a function of propagation distance (up to the 200th). The position of the centroid before and after the region of peak heating in each case follows a linear path, with deflection occurring predominantly within that region. By fitting a line to the trajectory before and after collapse, we can calculate the deflection of each pulse in terms of the deflection angle in milliradians (Fig. 8).

 

Fig. 7. Centroid location in the direction parallel to gravity of every 40th pulse, as a function of propagation distance.

Download Full Size | PDF

 

Fig. 8. Centroid deflection in milliradians as a function of pulse number (or time in ms).

Download Full Size | PDF

Direct comparison of our results to experiment is difficult, as the authors of [5] only include deflection results for Xenon at 2.7 atm. In simulation, the deflection in air reaches a steady-state value of less than −0.6 mrad. By comparison, for Xe at 2.7 atm the experimentally measured deflection reached approximately −3.5 mrad, with an apparent noise level of about $\pm ~ 0.5$ mrad. For the experiment performed in air, filamentation occurred in free space, as opposed to a closed gas cell, such that any deflection in this case was overwhelmed by random fluctuations [27]. One result which does lend itself to rough comparison, however, is Fig. 2(b) of [5], which shows the interferometrically-determined air density perturbation after a large number of pulses, 100 $\mathrm {\mu }$s before the next pulse. Figure 9 shows the simulated equivalent, after the 500th pulse. The simulated result compares qualitatively well with the experimental figure, which had a maximum density depression of $\delta \rho /\rho _0 \approx 0.065$.

 

Fig. 9. Simulated interferometrically-determined relative density perturbation $\delta \rho /\rho _0$, in analogy with Fig. 2(b) of [5]. This figure was calculated after 500 pulses, 100 $\mathrm {\mu }$s before the next pulse, using an interaction length of 9 cm.

Download Full Size | PDF

6. Conclusions

We have presented a theoretical model to describe the heating of air for ultrashort laser pulses and the interpulse thermal interactions of a pulse train in which the pulse separation allows for isobaric description of the hydrodynamics of air. The model includes the loss of laser energy to the air by means of linear absorption by molecules, and the nonlinear processes of photoionization, electron energization, and rotational Raman excitation. As such, the model can describe energy loss to the air for pulse durations from tens of femtoseconds to several picoseconds.

We calculated the air heating of a single laser pulse whose duration is much shorter than the acoustic transit time, and the resulting refractive index perturbation at isobaric time scales that would be encountered by subsequent pulses. We then calculated the cumulative index change over a laser pulse train where the primary cooling mechanism was thermal conduction.

We developed a 3D simulation in NRL’s PyCAP code based on this model to simulate the propagation of pulse trains undergoing nonlinear self-focusing. We compared the simulated peak electron density of a single pulse undergoing filamentation with a previous experimental study, and found that the simulation shows good agreement with experiment. We then self-consistently modeled the thermal interaction of 500 pulses undergoing filamentation in air at a repetition rate of 1 kHz. Later pulses were deflected downward by less than 0.6 milliradians, a result which could not be directly compared with experiment due to experimental noise issues but which is not inconsistent with the experimentally-observed level of noise. However, we compared our simulations with the interferometrically-measured air density perturbation and found agreement with experiment to within 40%.