Model for Ultrashort Laser-pulse Induced Ionization Dynamics in Transparent Solids

A comprehensive model of ultrafast laser-induced plasma generation intended for coupling with pulse propagation simulations in transparent solids is introduced. It simultaneously accounts for the changing spectrum of a propagating ultrashort laser pulse while coupling to the evolution of the energy-resolved nonequilibrium free-carrier distribution. The presented results indicate that strong pulse chirps lead to ionization dynamics that are not captured by the standard monochromatic treatment of laser-induced plasma formation. These results have strong implications for ultrafast laser-solid applications that depend on ionization in a strong nonlinear focus.


I. INTRODUCTION
Ultrafast nonlinear pulse propagation and laser-induced ionization in bulk media are interrelated research areas with broad and promising applications.A correct understanding of propagation effects combined with laser-material interaction is critical for progress in bulk micromachining [1,2], remote sensing [3], laser-induced breakdown spectroscopy, and laserbased medical procedures [4].To elucidate the physics of these applications, one must simultaneously model the propagation of the laser field and its interaction with the material [5].These combined areas of research present a challenge for theoretical calculations due in large part to the computing requirements of combining fully three-dimensional (3D) pulse propagation simulations with detailed calculations of laserplasma dynamics.Each of those calculations is cumbersome individually, and combining them typically involves using a detailed model of one process and a simplified model of the other.Such approaches therefore contain many inconsistencies and oversimplifications [6].
For most pulse propagation simulations in solids, the chosen models of photoionization and laser-plasma interactions were derived in the decades preceding ultrafast laser systems.As such, they are built on the assumption of plane-wave-like monochromatic radiation.The current approach is to replace the constant electric field amplitude of the monochromatic approximation with a time-dependent pulse amplitude, keeping the single-frequency approximation for the pulse spectrum [5,7].While this works well for pulses with narrow spectra, short pulses with appreciable bandwidth cannot be accurately represented in this way.Such an approach is also inconsistent with computational propagation schemes designed to fully account for arbitrary pulse spectra.
Laser-plasma models currently coupled to pulse propagation in solids also neglect many of the most important behaviors of carrier dynamics.This is particularly significant when the laser field drives the electron distribution far from its equilibrium configuration [8].Simulations for dielectrics frequently limit their laser-plasma models to a single rate equation for the conduction band electron density [7,9] or in some recent cases multiple rate equations for the same purpose [10][11][12][13][14].These models use phenomenological approaches for laser energy absorption by the carriers and are coupled to pulse propagation via a spatially local Drude model.Researchers specifically concerned with the evolution of these nonequilibrium distributions use quantum kinetic models to include many-body effects involving carrier collisions with laser photons, impurities, phonons, and other electrons [8,15,16].Currently, there are no such models of sufficient simplicity to couple with long distance (3+1)D pulse propagation calculations in solids.
To achieve greater fidelity to the fundamental physics of carrier scattering dynamics, pulse propagation simulations could model laser-plasma interactions in dielectric solids with models based on the quantum kinetic approaches often used to study semiconductors.This paper introduces such a model based on a set of free-carrier energy-resolved descriptions referred to as extended multirate equations.It is shown how to couple this laser-plasma model with pulse propagation and use the instantaneous pulse frequency to calculate chromatically representative ionization rates.Simulations using this approach reveal that pulse chirps that occur naturally during nonlinear propagation alter the nonequilibrium carrier distributions, thereby changing the optical properties of the plasma.Since previous simulations have shown that such effects can alter the structure of generated plasma channels when focused in a fused silica sample [17,18], the model therefore has strong implications for laser-induced ionization applications.Limitations of, and improvements to, the framework of the extended multirate equation approach are then discussed.

II. THEORY
In dielectrics and semiconductors, the free-carrier plasma results from electron transitions from the valence band to the conduction band, initially by photoionization.In bulk solids, this is typically modeled with monochromatic approaches, such as those derived by Keldysh [19] and others [20].Although widely used for ultrashort pulses, such monochromatic approaches do not properly represent the photoionization of multichromatic fields, such as one finds in strongly chirped pulses.Once in the conduction band, electrons may absorb laser energy and can transfer it to valence band electrons via collisions producing impact ionization.Sequential singlephoton absorption events result in a delayed development of the impact ionization as demonstrated by simulations using the Boltzmann [8] and Fokker-Planck [21,22] equations and then incorporated into single and multirate equations [10,12].These models of free-carrier dynamics to date are also monochromatic.
In most pulse propagation simulations, the time-evolved carrier density, calculated from a rate-equation model [9], is exported to a Drude model for the free current density that assumes averaged values for material constants, such as collision rates and effective carrier masses [7,[23][24][25][26][27].Since collision rates and effective masses are known to vary significantly as functions of electron energy [28], the free-carrier transport must sensitively depend on both the spectral components of the laser pulse and the conduction band energy (or momentum) at which laser photons are absorbed or emitted.Furthermore, both experiment and pulse propagation simulations demonstrate that the energy absorbed by the plasma from ultrashort pulses varies as a function of the instantaneous intensity slope and pulse chirp [17,29,30].A self-consistent theoretical framework to address these issues is introduced in the following.All expressions are in SI units.

A. Photoionization
In this paper, the photoionization rate derived by Keldysh for solids is used [19].The Keldysh formula for a complex field E with angular frequency ω in a solid of band gap U and reduced electron-hole mass m r is [19] Here, the Keldysh parameter and E 2 (x) are complete elliptical integrals of the first and second kind, respectively, as defined in Ref. [31].The function Q(γ,x) is given by . .denotes the integer part, and (z) = z 0 exp(y 2 − z 2 )dy is the Dawson function.The relationship of the Keldysh formula to optical intensity and frequency is shown graphically in Fig. 1.One can view this plot as comprising many ridges, each of which corresponds to a range of frequencies and intensities where a particular order of multiphoton ionization occurs.In other words, a region where W PI ≈ σ k (ω)I k where σ k is the kth-order multiphoton ionization coefficient and I is the optical intensity.For example, the ridge present between 300 and 400 nm at low intensity is a three-photon absorption region, whereas the ridge to the right of it (beginning just after 400 nm at low intensity) is a four-photon absorption region, etc. Interesting laser-plasma dynamics can be expected for cases where the spatiotemporal intensity and/or spectral content of a given laser pulse spills over from one multiphoton ionization order to another.Figure 1 demonstrates the strong frequency dependence in the visible to near-IR range.However, note that at much longer wavelengths (in the tunneling regime not shown in the plot), this frequency dependence subsides [32].
When calculating the photoionization rate, the instantaneous frequency ω(t) = −Im[(∂ t E)/E] is used in Eq. ( 1), where the notation Im[. ..] denotes the imaginary part and E(t) is the complex forward-propagating field.Using the instantaneous frequency in the Keldysh formula allows one to approximate the effect of a changing pulse frequency, which occurs naturally during nonlinear propagation due to self-phase modulation and other optical processes.For a negatively chirped laser pulse, the blue-shifted photons on the leading pulse edge can lead to higher photoionization rates earlier in the laser pulse.This in turn allows energy in the trailing pulse edge to contribute more strongly to the total ionization yield through impact ionization.The opposite situation can occur for a positively chirped laser pulse.The reality of this effect was demonstrated experimentally when it was shown that the surface damage threshold for fused silica was up to 20% lower for a negatively chirped ultrashort laser pulse, as compared to that of an otherwise identical positively chirped pulse [30].

B. Ionization dymanics
Once free carriers are generated, they may interact with the surrounding ions, the external laser field, and with each other.To model these dynamic processes, this work proposes using a coupled system of extended multirate equations (EMREs) that discretize the carrier distribution and free-current density in energy-space [see Eqs. ( 2) and (16), respectively].Multirate equations (MREs) for the carrier distribution have been proposed in the past [11,12] and successfully used in pulse propagation simulations [24].The MRE, in its originally proposed form [12], discretizes the carrier distribution into conduction band energy regions of = ω 0 , where ω 0 is the frequency of the applied laser field.In the MRE model, it is assumed that the laser field frequency does not change appreciably during exposure.The MRE approach was later extended for use in semiconductors by further discretizing the energy space to allow energy deposition from the free-carrier plasma into the lattice [10], hence an "extended" multirate equation.
The proposed EMRE for dielectrics is expressed as where i is the ith represented conduction band energy and n i is the density of free carriers in the energy range i + = i+1 .The energy spacing should be equal to or smaller than the smallest energy transition the calculation is to resolve.The terms on the right-hand side of Eq. (2) describe contributions from photoionization, impact ionization, electron-phonon collisions, electron-electron collisions, and one-photon absorption/emission by free carriers, respectively.As with any MRE, the total number of carriers per volume is given by N = i n i and the total energy per volume of the plasma by E N = i n i i .
In general, all terms on the right-hand side of Eq. ( 2) should be calculated in the manner that best represents the chosen material.A first-principles-based quantum mechanical approach for all terms is naturally preferable [32].However, when coupling to multidimensional pulse propagation simulations this may not always be possible in practice.In anticipation of such cases, this paper assumes simplified descriptions that capture the underlying variation of electron-energy-dependent processes, but for which the computation effort required is comparatively small.The photoionization term is taken to be where is the energy at which photoionized electrons enter the conduction band, ˜ is the effective band gap [19], and ω(t) is the instantaneous frequency.The contribution of impact ionization is given by where is the critical energy for impact ionization, μ cv is the reduced effective mass of the conduction and valence electrons, g = ˜ / , m = max / , and max is the highest modeled conduction band energy [10].The avalanche contribution allows electrons with energy exceeding c to impact a valence electron and promote it to the conduction band, losing an energy equal to the effective band gap ˜ in the process.The energy dependence of the impact ionization rate α( i ) is approximated by the so-called Keldysh impact ionization formula , where P imp is the impact rate coefficient [8] and θ (x) is the Heaviside step function.
The electron-phonon term can be approximated by a net plasma energy relaxation into the phonon gas and is given by ṅe-pn Here, τ pn i is the energy-dependent electron-phonon scattering time and the index p = pn / .The scattering time is implemented in terms of the scattering rate ν pn i and the limits of the energy discretization: , where pn is the mean phonon energy.The total energy deposited from the plasma into the lattice can then be calculated by pn i (n i /τ pn i )dt, from which one can calculate the temperature of the lattice.
The electron-electron collisions are also represented by a relaxation approach ṅe-e i = Here, n eq i (t) is the energy-resolved equilibrium distribution to which the carriers would relax provided only electronelectron collisions were allowed.For an insulator with low conduction electron densities or high temperature, the Fermi-Dirac statistics of conduction electrons may be approximated by a Maxwellian [33,34], and thus the resulting equilibrium distribution is given by Here, g( is the density of states per volume for a free-electron gas, T is the temperature of the equilibrium distribution, and μ c is its conduction band referenced chemical potential.The temperature and chemical potential of n eq ( ,t) are determined at each time step prior to constructing the distribution.This is done in such a way as to ensure that the total number of particles per volume is conserved ( i n i = i n eq i ) as well as the total energy per volume ( i i n i = i i n eq i ) during elastic electron-electron collisions.Following the approach of Ref. [35] the Debye screened electron-electron collision rate is calculated as where avg = E N /N is the average energy per carrier.The last term on the right-hand side of Eq. ( 2) is given by The terms W abs i−j and W ems i+j represent the carriers per volume per time entering the i region by one-photon absorption and emission.Here, carriers entering an energy region represented by i do so by absorption or emission of field energy from carriers at energies i+j = i + ω(t) or i−j = i − ω(t), respectively.Correspondingly, the rates W abs i and W ems i give the total carriers per volume leaving the i region by onephoton absorption and emission events.From basic EM theory, the work done per volume per time on the charges by the laser field is J i (t) • E(t).In this paper, the time average of this quantity is used to determine the one-photon absorption/emission contributions to Eq. ( 2): where ω(t) is again the instantaneous frequency.
When W i is positive, the result is interpreted as the net rate of carriers per volume promoted to an energy ω(t) above i .However, it is important to note that the quantity W i can be negative.A negative value of W i can be understood as the plasma depositing a net energy back into the field (emission).Thus, the absorption and emission rates are calculated by There are two primary advantages of this EMRE compared to an ordinary MRE (or a single rate equation) for the carrier density when coupling to ultrashort pulse propagation.(i) It accounts for the multichromatic character of an ultrashort pulse spectrum undergoing nonlinear propagation.The use of the instantaneous frequency provides a mechanism for calculation of photoionization and one-photon absorption in any optical setting where the instantaneous frequency is meaningful.(ii) It provides a modular framework for investigating the effects of individual solid-state processes dynamically coupled to multichromatic sources.If a specific solid-state process is of particular interest, then one of the models described above can be replaced with an alternative, more sophisticated model without affecting the general framework of the EMRE.

C. Free-current density
Rate-equation models for conduction plasma densities to date have calculated energy absorption from the laser field under the monochromatic approximation.When coupled with pulse propagation simulations, the resulting carrier density is then used to calculate a free-current density for use in the propagation equations, typically by a Drude model [3,7,36].One consequence of this approach is that laserenergy absorption appearing in the pulse propagation model is not necessarily consistent with the laser-energy absorption appearing in the plasma evolution model.A primary goal of the proposed approach is to force agreement between the laser-material model and the pulse propagation model by a coupled calculation of the one quantity they both share in principle; the free-current density.
The evolution of the free-current density must be directly coupled to Eq. ( 2), and therefore identically discretized in energy space.The total current density can be expressed as J = qN v avg , where q is the free-carrier charge, N is the number of carriers per volume, and v avg is the average velocity of the plasma "fluid."What is desired is a current contribution J i from each population of carriers at average energy i such that the total free-current density is given by A current density can also be expressed in terms of the population momentum per volume of the carriers: J i = (q/m i ) p i , where p i = m i n i v i , m i is the effective mass, and v i is the average velocity of the carriers at energy i .For an isotropic material, the net momentum of any carrier population is assumed to be zero in the absence of an external field.Once exposed to a field, the population begins to oscillate, gaining and losing momentum from the field.Therefore, while Eq. ( 2) can be understood by tracking changes in carrier number and energy, an EMRE for the current density should be derived by tracking the corresponding changes to the carrier momentum.For this purpose, it is useful to reorganize the terms in Eq. (2) as follows: Here, ṅ0 i contributions add to the carrier population at energy i but do not add any average momentum to those carriers (i.e., they enter the population at zero average velocity [3]), ṅ− i contributions reduce the carrier population at energy i and take their momentum with them, and ṅ+ i contributions add to the carrier population at energy i and bring momentum with them.
To model the change in total carrier momentum due to electron-phonon-photon scattering, this work uses an energydependent momentum relaxation time of τ c i .This contribution represents the only loss mechanism included in this paper.However, note that one may easily include others such as those from scattering of electrons with ions or electrons with neutral atoms [35], assuming that these effects are simultaneously included in Eq. (2).Using this approach, the momentum relaxation time provides a "causal" numerical mechanism from which the laser energy absorption term in Eq. (2) will be calculated.This approach, while self-consistent numerically, is admittedly backwards physically.Physically speaking, it is the absorption/emission terms in Eq. ( 2) that should give rise to momentum loss in the total carrier population.If one calculates the laser absorption/emission terms in Eq. ( 2) from purely quantum methods, then those terms will provide the causal numerical mechanism for momentum loss for a current density EMRE, instead of the approach taken here.
Combined with the change in momentum due to the external laser field, a corresponding set of equations representing the total change in momentum per volume of each population is given by Here, the total momentum per volume of the free carriers is p = i m i n i v i and the average carrier momentum is defined as p avg = p/N.It is assumed that populations leaving the energy i by electron-electron collisions take their momentum with them, while carriers entering that level by electronelectron collisions enter with the average carrier momentum.Using these expressions, the time rate of change for the current density must be where the effective rate 1/τ eff i is defined as The current J ee i = qn eq i p avg /m i is the contribution of incoming carriers from electron-electron collisions.Equation ( 16) allows the calculation of absorption and optical effects from carriers of every energy region.Once the solutions to this EMRE are summed [as per Eq. ( 13)], one has the free-current density derived directly from the nonequilibrium distribution including multichromatic (pulse-chirp) effects.
Note that in the case of an approximately constant effective mass m i ≈ m const , summing Eq. ( 16) over all energies gives the form of a traditional "Drude" model: The last term in Eq. ( 18) gives the damping contribution to the total current.Note, however, that the current density cannot be solved in the form of Eq. ( 18) without significant assumptions about the nonequilibrium distribution.Rather, Eq. ( 16) should be used in order to capture the full energy dependence of the electron distribution on the damping term.

III. SIMULATIONS
In this section, the implementation of Eqs. ( 2) and ( 16) is discussed and results for three ultrashort pulses of different initial chirps are presented.Material parameters are chosen to correspond to fused silica and are summarized in Table I.
Also common to all simulations in this paper are the electron-phonon energy and momentum relaxation rates as a function of electron energy.While these would ideally be calculated using a quantum mechanical or semiclassical method, the intention for the proposed EMREs is to couple them to long-distance propagation models.Therefore, a more pragmatic approach is taken here.Following the example of Ref. [9], this work extracts the energy and momentum relaxation rates from Ref. [28] that were calculated for fused silica by solving a quantum mechanical Boltzmann transport model.The corresponding relaxation times are plotted in Fig. 2. Note that the electron energy relaxation rate calculated in Ref. [28] contains contributions from impact ionization which must be subtracted away before producing the electronphonon time in Fig. 2.
The pulse parameters for all simulations reported in this paper are summarized in Table II and are chosen to correspond approximately to a 60-fs frequency-doubled Ti:sapphire laser field.The laser pulse is linearly polarized and is represented by a scalar field E(t) for which the spatial dependence is neglected in this study.The field is assumed to be Gaussian in time and is numerically constructed according to the following expression: where ω 0 = 2πc/λ is the central frequency of the pulse.The chirp parameter β takes on either a zero value for the unchirped case or a positive (negative) value for a positive (negative) linear chirp.For the chirped cases (β = ±2 × 10 27 fs −2 ), note that this combination of pulse width and linear chirp is equivalent to that of 8.1-fs transform-limited pulse stretched to 60 fs.During nonlinear propagation of ultrashort pulses, strong positive chirps of this kind develop naturally due to self-phase modulation, assuming the nonlinear refractive index to be positive.The unchirped pulse (β = 0) is transform limited at 60 fs and, although identical in intensity and pulse width to the chirped pulses, is spectrally district from the other two pulses.Note also that the standard methods of calculating plasma generation for propagation simulations make no distinction  between these pulses when calculating ionization yields.They instead assume that such calculations depend only on the instantaneous intensity and central frequency, thus obtaining the same unchirped result for all three cases.The remainder of this section explores the resulting ionization yields and optical interactions generated by the two pulses of opposite chirps.The additional inclusion of the nonchirped pulse case further demonstrates how neglecting multichromatic effects influences the calculation.

A. Role of photoionization
The calculated photoionization rates for each pulse are shown in Fig. 3.Note that the rate for the unchirped pulse, shown as the solid black line, is the monochromatic rate that would be commonly used for each of the pulses in pulse propagation calculations.The fact that each chirped pulse leads to a peak photoionization rate that is greater than the peak rate of the unchirped pulse will be important.This will provide a higher number of "seed" electrons for the avalanche process.However, it is more important when those individual peaks occur in time.For the case of the negatively chirped pulse, the peak rate occurs on the leading pulse edge.This leads to more avalanching during the duration of the pulse, maximizing the ionization yield.For the unchirped and positively chirped pulses, it is not obvious which case will yield more ionization events.The unchirped pulse has a lower overall peak photoionization rate, but its first peak rate occurs earlier in the field history.The relative efficiency of avalanching therefore is decisive for determining which case (positive chirp or no chirp) will cause greater total ionization.

B. Solutions to the EMREs
Figure 4 shows the solutions to Eq. ( 2) as a function of electron energy and time for the three variously chirped pulses described in Sec.III.Note that Fig. 4 plots a more traditional quantity for representing the electron distribution; the number of electrons per volume per energy, f ( ,t) = dn/d .This quantity is obtained from solutions to the EMRE by dividing by the energy spacing .The corresponding total ionization yields and average carrier energies as functions of time for each pulse are shown in Fig. 5.The pulse is centered about t = 0. Note that the nonequilibrium behavior of the distribution for each pulse appears during the initial ionization events before t = 0.After the peak of the pulse has passed, all the distributions rapidly move towards a quasiequilibrium on the time scale of tens of femtoseconds.Here, the term "quasiequilibrium" means only that the plasma has relaxed to a smooth, continuous, and slowly evolving distribution, but is not yet in thermal equilibrium with the background lattice.This result is in good agreement with the theoretical results of Kaiser [8] which solved an energy-resolved quantum Boltzmann equation for electrons and phonons in fused silica exposed to ultrafast laser fields.One major factor determining the nonequilibrium behavior in Fig. 4 is the initial energy of photoionized electrons, given by Eq. ( 4).This initial placement of electrons, which changes as the intensity and instantaneous frequency changes, determines how effectively the carriers initially absorb and scatter the laser field according to the momentum relaxation times in Fig. 2. Through this effect, the instantaneous frequency in the proposed model significantly impacts not only the photoionization rate, but also the optical absorption and scattering generally associated with a Drude model.Using Eq. ( 18), one can obtain a distribution-averaged momentum relaxation time as a function time for each pulse.This result is shown in Fig. 6(a) for the time window where the field intensity is appreciable.
Figure 6(a) demonstrates how multichromatic pulse effects alter distribution-averaged relaxation time, often assumed to be a constant of the material.These multichromatic effects occur primarily for the leading edges of the pulse in time because the distribution has a strong frequency-dependent nonequilibrium configuration there (see Fig. 4).On the trailing pulse edge the distributions are all converging on quasiequilibrium configurations with the same average energies despite having distinct ionization yields spanning an order of magnitude (see Fig. 5).Here, all three cases approach a value on the order of 0.4 fs, which is in good agreement with the experimentally fitted value of this parameter by Gigúere et al. [37].
In the results presented here, relaxation to the quasiequilibrium is primarily due to electron-electron collisions.The calculated electron-electron collision times, defined by Eq. ( 9), are plotted for each pulse in Fig. 6(b) as functions of time.It is notable from Fig. 6 how similar the forms of τ c and τ ee actually are.This is somewhat surprising given that, in the proposed model, τ ee was not used directly in the calculation of τ c i , although the authors note that this is commonly done as per Matthiessen's rule [10,35].This is simply an indication of each parameter's dependence on the average energy avg of the distribution.Also, this result does not support the occasionally assumed form of a carrier-density-dependent Drude collision time, most commonly τ c ∝ N −1 , since the τ c for all pulses converges to a single value despite N being different by an order of magnitude.In particular, it is notable that the averaged momentum relaxation times in Fig. 6(a) are all well approximated by the formula τ c ≈ (1.4 eV fs)/ avg .

C. Pulse propagation
To investigate the effects of the initial chirp on pulse evolution, 1D propagation simulations are performed for a 20-μm penetration depth in fused silica.These simulations solve Eqs. ( 2) and ( 16) simultaneously with a unidirectional pulse propagation equation (UPPE) for the forwardpropagating electric field in the frequency domain as it evolves along propagation axis z [6]: Here, is the Fourier transform frequency coordinate, k( ) = n( )/c is the frequency-dependent wave vector, and n( ) is the linear index of refraction as calculated by a Sellmeier equation.Equation ( 20) is a fully spectral ODE accounting for linear and nonlinear dispersive effects to infinite order [7].Suppressing the z dependence for notational simplicity, the nonlinear polarization P NL ( ) contribution is calculated in the time domain by [6]  where n 0 = n(ω 0 ), 0 is the permittivity of free space, n 2 is the nonlinear refractive index, and I (t) is the intensity.The nonlinear response function R(t) is given by [29] e −t/τ 2 sin(t/τ 1 )θ (t), (22) where f r is the fraction of the nonlinear response attributable to the nuclear motion, τ 1 and τ 2 are the Raman response times, and θ (t) is the step function.The values of the nonlinear propagation constants used are n 2 = 3.2 × 10 −16 cm 2 /W [38], τ 1 = 12.2 fs, τ 2 = 32 fs, and f r = 0.18 [24].The current-density term in Eq. ( 20) is also calculated in the time domain.It comprises a photoionization current J PI (t) as well as the free-current density J f (t) [7,29]: The free-current density is provided by solving Eq. ( 16).The second current resulting from the photoionization absorption is given by [39,40] where l PI = ˜ / ω(t) + 1 [8].Note also that alternative treatments of the ionization current exist and, at intensities of 10 14 W/cm 2 and higher, can lead to a dephasing from the free-current density [41].
The simulations are performed for the three chirped pulses described in Sec.III.The incident field is numerically constructed according to Eq. (19).Upon interacting with the sample surface, a field transmission coefficient of T = 2/[1 + n eff ( )] is applied to the field in the frequency domain where the effective index n eff ( ) is calculated by n 2 eff ( ) = 1 + P eff ( )/ 0 E( ), where the effective polarization is defined as P eff ( ) = 0 [n 2 ( ) − 1]E( ) + P NL ( ) − J ( )/i .Figure 7 shows the peak ionization yield and transmitted pulse energy as a function of propagation distance in the medium.Since the simulations are 1D, Fig. 7(b) technically shows the transmitted fluence, the initial fluence having a peak value of 1.6 J/cm 2 .As first seen in Fig. 5, the initial ionization yield resulting from the negatively chirped pulse is nearly an order of magnitude higher than those of the positively chirped and unchirped pulses of identical pulse width and energy.The negatively chirped pulse continues to produce higher ionization yields than the other pulses for about 20 μm.This has a significant effect on the total absorption of laser energy as seen in Fig. 7(b).Excluding the 4%-5% of energy that is reflected at the surface, the unchirped and positively chirped pulses result in about 20% absorption whereas the negatively chirped pulse results in about 45% absorption.Thus, the negative pulse chirp doubles the total energy absorption in this case.
In Fig. 8, the instantaneous frequency, intensity, and spectra of the three pulses are shown.The left column of this figure shows plots for the incident pulses while the right column shows the same pulses after 20 μm of propagation.Differences between the incident and final instantaneous frequencies show the standard effect of self-phase modulation, which imparts an intensity-dependent positive chirp for positive values of n 2 .The differences between incident and final pulse intensities show the influence of ionization on the pulse shape, which tends to absorb the trailing edge while leaving the leading edge intact.Comparison of the initial and final spectra reveal that plasma-induced blue-shifting occurs in all cases, as shown from earlier works [40][41][42].However, for the positively chirped pulse, there is more red-shifting than blue-shifting of the pulse energy.This occurs because the pulse is red-shifted on the leading edge [see Figs.8(a) and 8(b)], which is not as strongly affected by free-carrier absorption.This fact combined with the additional positive chirp from self-phase modulation causes a net red-shift of the pulse spectrum.The opposite occurs for the negatively chirped pulse, which is the most strongly blue-shifted of all as a result.It should be noted that the 1D simulations presented here make no account of self-focusing or plasma defocusing effects.When the transverse spatial dependence of a finite laser beam is taken into account by higher-dimensional simulations, additional complexity arises in the propagation dynamics through behaviors such as filamentation [3].

V. DISCUSSION
The presented results approximate the influence of pulse chirps on laser-induced ionization and free-carrier optical effects.The pulses were chosen to have the same pulse widths and total energy but different time-frequency arrangements.The multichromatic effects of such pulses enter model primarily through the influence of photoionization.They do this in two key ways: (i) by changing the magnitude and temporal position of the peak photoionization rate, thus providing more (or less) time for avalanching to increase the ionization yield and (ii) by virtue of their initial placement on the conduction band energy scale, this placement varying as a function of intensity and instantaneous frequency.Combined with the conduction band energy-dependent description of carrier relaxation rates, the subsequent evolution of the plasma's optical properties is strongly influenced by these initial carrier distributions.These effects on the optical properties are most strongly seen on the leading edge of the pulse.On the trailing edge of the pulse, the multibody collisions, electron-electron interactions in particular, move the distribution rapidly towards quasiequilibrium on a time scale of tens of fs or less.Regardless of initial pulse chirp or the total ionization yield, all the resultant distributions in Fig. 4 evolved toward the same average energy and thus nearly identical relaxation times on the trailing edge (see Figs. 5 and 6).
There are some notable limitations to the presented approach as well as opportunities for future improvement.First, the proposed approach relies on the instantaneous frequency both in calculating the photoionization rate as well as conduction band free-carrier absorption.The Keldysh photoionization rate was not derived with such an assumption, anymore than it was derived for fields of changing intensity, for which it is ubiquitously used.Also, the instantaneous frequency is not always a meaningful or representative quantity, particularly in the cases where the temporal phase is difficult to calculate accurately.For example, this will occur in pulse propagation simulations when events such as pulse splitting cause small, but numerically noisy, field amplitudes between the two split pulses as they interfere.In such situations, it would be prudent to use a different "average frequency" value.In any event, photoionization in these small-amplitude time regions is comparatively small.Alternatively, for photoionization there have been recent attempts to treat pulses of arbitrary temporal shape and phase [43][44][45].The authors note that the method derived in Ref. [44] is particularly well suited to an EMRE framework.
Another limitation of the presented approach is its reliance on pretabulated relaxation rates for a particular material.A possible improvement is to use self-consistent quantum approaches to calculate the photon absorption rates and electronphonon collision rates for the carriers [46].Furthermore, the influence of many-body collisions on carrier distributions is manifest in Fig. 4.These results also rely heavily on the relaxation approximation.In particular, the assumed models of impact ionization and electron-electron collisions have a strong influence on the shape of the carrier distribution, and hence the energy-averaged optical properties.These processes would be more accurately treated with a statistical approach.This would enforce a more physically realistic relaxation to quasiequilibrium from electron-electron collisions.For impact ionization, this would allow impacted and impacting electrons to occupy all energies subject to the constraint of conserving both total momentum and energy [8].

VI. CONCLUSIONS
In summary, a comprehensive model of ultrafast laserinduced plasma dynamics intended for coupling with pulse propagation simulations in transparent solids was introduced.The model comprised two sets of extended multirate equations for the free-carrier distribution and the free-current density.The changing spectrum of a propagating ultrashort laser pulse is taken into account by use of the instantaneous frequency in calculating photoionization and single-photon absorption rates.It simultaneously couples to the evolution of the energyresolved nonequilibrium free-carrier distribution.Simulations solving this model indicate that strong pulse chirps, such as those occurring naturally during nonlinear propagation, lead to ionization dynamics that are not captured by the standard monochromatic treatment of laser-induced plasma formation.The order-of-magnitude chirp-dependent differences in ionization yield and the electron distribution shapes govern the strength and character of the laser-plasma interactions.The proposed model therefore provides a framework for laser-material interactions that is practical and provides better insight into the physics of laser-induced ionization by intense ultrashort laser pulses.

FIG. 1 .
FIG. 1. (Color online) Keldysh's photoionization rate as a function of intensity and laser wavelength in fused silica with a band gap of 9 eV and reduced electron-hole mass of 0.86 rest electron masses.

FIG. 3 .
FIG. 3. (Color online)Photoionization rates as a function of time for the unchirped, negatively chirped, and positively chirped pulses according to Eq. (1) using the instantaneous intensity and instantaneous frequency.

FIG. 4 .
FIG. 4. (Color online) Electron energy distributions after exposure to (a) the unchirped pulse, (b) the negatively chirped pulse, and (c) the positively chirped pulse described in Sec.III.

FIG. 5 .
FIG. 5. (Color online) (a) Total ionization yield N and (b) the average electron energy of the distributions from Fig. 4 as functions of time.

FIG. 6 .
FIG. 6. (Color online) (a) The average momentum relaxation time τ c and (b) the electron-electron collision time τ ee of the distributions from Fig. 4 as functions of time.

FIG. 7 .
FIG. 7. (Color online) (a) The peak ionization yield N max and (b) energy transmission (pulse energy in units of the incident pulse energy) as functions of the propagation distance z.

FIG. 8 .
FIG. 8. (Color online) The instantaneous frequency (a) and (b), intensity (c) and (d), and normalized transmitted spectra (e) and (f), at propagation distances of z = 0 and 20 μm, respectively.The negatively chirped pulse is shown as the blue dashed line, the unchirped pulse is shown as the black solid line, and the positively pulse is shown as the red dotted line.

TABLE I .
Material parameters used in simulations.

TABLE II .
Pulse parameters used in simulations.