Numerical simulations of ultrashort laser pulse ablation and plasma expansion in ambient air

Spectrochimica Acta Part B 56 Ž2001. 973᎐986

Numerical simulations of ultrashort laser pulse ablation and plasma expansion in ambient air 夽 b , M. Chaker a , F. Vidal a,U , S. Laville a , T.W. Johnston a , O. Barthelemy ´ B. Le Drogoff a , J. Margot b, M. Sabsabi c

´ Institut National de la Recherche Scientifique-Energie et Materiaux, 1650 Boul. Lionel-Boulet, Varennes, Canada, ´ QC J3X 1S2 b Departement de Physique, Uni¨ ersite´ de Montreal, ´ ´ Montreal, ´ C.P. 6128, Succ. Centre-Ville, Montreal, ´ Canada, QC H3C 3J7 c Institut des Materiaux Industriels, Conseil National de Recherches Scientifique du Canada, 75 Boul. de Mortagne, ´ Boucher¨ ille, Canada, QC J4B 6Y4 a

Received 8 November 2000; accepted 21 March 2001

Abstract Using a self-consistent one-dimensional Cartesian Lagrangian fluid code, we modeled the ultrashort laser pulse ablation of solid aluminum and the subsequent plasma expansion in ambient air. A laser fluence of approximately 10 Jrcm2 is considered. The code axial plasma temperature and density are strongly inhomogeneous and the maximum radiation emission generally occurs in the front of the plasma. The code average plasma temperature is in good agreement with the experiments for all times, while larger discrepancies with respect to the experiments are observed at late times for the plasma density. Experimental results are in reasonable agreement with the condition of thermodynamic equilibrium, which is an important assumption in the model. 䊚 2001 Elsevier Science B.V. All rights reserved. Keywords: Laser; Ablation; Plume; Modeling; Plasma

夽

This paper was presented at the 1st International Congress on Laser Induced Plasma Spectroscopy and Applications, Pisa, Italy, October 2000, and is published in the Special Issue of Spectrochimica Acta Part B, dedicated to that conference. U Corresponding author. E-mail address: [email protected] ŽF. Vidal.. 0584-8547r01r$ - see front matter 䊚 2001 Elsevier Science B.V. All rights reserved. PII: S 0 5 8 4 - 8 5 4 7 Ž 0 1 . 0 0 1 9 5 - 1

974

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

1. Introduction Laser induced plasma spectroscopy ŽLIPS. is used successfully in the determination of the relative concentrations of the components of various solid samples, including metallic alloys, insulators, archeological artifacts, and paintings Žsee w1,2x for a review of the LIPS applications.. There is, however, a constant demand for improving the precision of the method and, in particular, for the detection of the smallest amounts possible of chemical elements present in the sample to analyze. This translates, in practice, in obtaining the highest signal-to-noise ratio, i.e. the best resolution of the spectral lines over the background or continuum emission. For this purpose, it is natural to investigate the influence of the various laser parameters Žpulse length, wavelength, fluence, etc.. and of the ambient gas conditions Žnature and pressure.. In this respect, our recent experiments suggest that ultrashort Ži.e. subpicosecond. laser pulses can achieve a better line resolution at early times than longer pulses w3x. With the aim of improving our understanding of the physical processes involved in LIPS and of the role of the various control parameters, we have undertaken Žas well as an experimental effort w3x. the development of a computer code treating both the processes of laser ablation and subsequent plasma expansion in a given ambient gas. The code is a one-dimensional fluid model including a self-consistent treatment of hydrodynamics, laser energy absorption, thermal conduction, and realistic thermodynamic properties. A particular feature of the model used here is that there are no adjustable parameters. Similar codes have been used in the past few years to model the phenomena of ultrashort laser absorption w4᎐7x, shock waves propagation in solids w8᎐11x, and laser ablation of metals w12x and of semiconductors w13x. However, to our knowledge, such codes have not been used to simulate the properties of the expanding plasma in an ambient gas. Different codes have been developed to model the expansion of the ablated matter only w14᎐20x. In those works, however, the initial conditions were defined more or less arbitrarily. In this paper, we present simulation results for

ultrashort laser pulse ablation and the subsequent plasma expansion in ambient air, and compare them to our experimental results. Because its properties are so well known, aluminum is considered here. However, we expect our results to apply, at least qualitatively, to other conductors. Besides the mentioned-above fact that ultrashort laser pulses seem to be promising for improving the precision of the LIPS technique, they are fortunately easier to deal with from the point of view of modeling. For ultrashort laser pulses, the matter does not have time to expand within the pulse length and most of the laser energy is thus absorbed at solid density. In that case, the mechanisms of laser energy absorption are relatively well understood w4,5x. For longer pulses, more laser energy is absorbed in the low-density expanding matter, where the absorption mechanisms are often more complex w21x. Sooner or later, the presence of an ambient gas involves a further degree of complexity in the model due to the strong interactions between the gas and the expanding plasma. Since the initial plasma velocity is of the order of the sound velocity in the cold metal, a violent shock wave is generated in the ambient gas Žwhere the sound velocity is much smaller. which resists at first feebly, then more strongly, the plasma expansion and eventually stops it completely. In this work, we consider air at normal pressure as the ambient gas. This article has two objectives. The first one was to show, through comparisons with experiments, that our model provides a good description of the phenomena that are relevant for LIPS Žat least for the limited set of parameters considered here.. The second objective is to provide, through the examination of the code results, some insights of the expanding plasma properties as a function of time and space. A detailed investigation of the role of the pulse length Žand other laser parameters. and of the ambient gas conditions on the expanding plasma are beyond the scope of the present article and will be presented elsewhere. After having presented the model used in the next section, we present, in Section 3 the simula-

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

tion results for laser ablation and plasma expansion. The conclusion is presented in Section 4.

2. The model In our computer model, the fluid equations for the conservation of mass, momentum and energy are solved by means of a one-dimensional Lagrangian code Žfor the Cartesian coordinate along the normal to the solid surface, z . whose numerical scheme is essentially based on the code MEDUSA w22x. We treat the plasma as a single fluid characterized by one velocity u and one temperature T. In the Lagrangian scheme, the surface mass of each cell is constant in time. The velocity u is obtained by solving the Navier᎐Stokes equation: ␳

du ⭸ Ž pq q . sy dt ⭸z

Ž1.

where t is time, ␳ is the aluminum mass density, p is the total pressure and q is the von Neuman artificial viscosity. The temperature is obtained from the energy conservation equation:

ž ⭸U ⭸T /

␳

dT q dt

⭸U ⭸␳

žž /

T

y

p ␳2

/

d␳ sS dt

Ž2.

where U is the local internal energy. The source term is S s X q H q R q J where X is the rate of absorption of the laser radiation, R is the viscosity term, and H s y␳y1 Ž⭸Kr⭸z . is the net heat conduction contribution. The heat flux K itself is given by the model of Lee and More w23x, which holds for a large range of density and temperature, including the cold solid conditions and the hot tenuous plasma. For simplicity, we have neglected the rate of radiation loss J. The laser energy absorption X is treated by numerically solving the Helmholtz wave equation w24x for the complex amplitude of the electric field E in S-polarization, and incident angle ␪ with respect to the normal to the solid surface, d2 E q k 02 Ž ␧ y sin 2 ␪ . Es 0 d z2

Ž3.

975

where k 0 is the laser wavenumber in the vacuum and ␧ Ž ␻ . s 1 q i4 ␲ ␴ Ž ␻ .r␻ Žcgs units. is the complex dielectric function. In this expression, ␻ is the laser angular frequency and ␴ Ž ␻ . s ␻2pr4␲ Ž ␯ei y i␻ . is the complex electrical conductivity in the Drude approximation. The electronion collision frequency is ␯ei s ␻2pr4␲ ␴0 where ␻ p is the plasma frequency and ␴0 is the d.c. conductivity, also obtained from the model of Lee and More w23x in a way consistent with the thermal conduction. The rate of laser energy deposition is then X s ReŽ ␴ .< E < 2r2. We have checked that the results of our code are in good agreement with other work w4,5x for the laser energy absorption in aluminum. For the calculation of the pressure p and internal energy U of the solid aluminum and the subsequent plasma, we use the quotidian equation of state ŽQEOS. model w25x which combines the Helmholtz potential for free electrons ŽThomas᎐Fermi model., for ions ŽCowan model. and for the cohesion forces ŽBarnes model.. ŽThe latter are adjusted to cancel the pressure of the cold solid.. QEOS provides a well-behaved approximate equation of state for any material for a given material solid density, atomic number and bulk modulus. Despite the fact that QEOS can easily handle a separate treatment of ion and electron temperatures w10x, the required treatment of the electron-ion thermal coupling would here only add inessential complications. The shock wave in ambient air created by the expanding plasma can, in principle, be dealt with using the same formalism as for aluminum. Unfortunately, we have found that such an approach is very costly in computer time because an increasing number of meshes in the calculation grid is required as the shock propagates. Moreover, the time step must be kept small for the code to remain numerically stable. Consequently, we adopted a simplified approach, which consisted of treating the shock wave in air as if it were always in equilibrium with the expanding plasma. Furthermore, we have neglected the phenomena of molecule dissociation and of ionization in the shock wave, as if the shock was weak. Using the Hugoniot relations Žconnecting the upstream and

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

976

downstream regions of the shock wave. w26x, one obtains the shock pressure: ps s p 0 1 q

ž

2␥ Ž M 2 y 1. 1q␥

/

Ž4.

where p 0 is the ambient air pressure and ␥ is the usual specific heats ratio which is 9r7 in cold air. Defining Vp as the velocity of the plasma᎐air boundary, the Mach number is, for Vp ) 0: Ž 1 q ␥ . Vp Ms 1q 4 c0

ž žž

4 c0 Ž 1 q ␥ . Vp

1r2

2

/ / q1

/ Ž5.

while Ms 1 for Vp F 0. In this expression, c0 is the sound velocity in ambient air. The velocity Vp is obtained by solving the Navier᎐Stokes equation wEq. Ž1.x, where the shock pressure ps is used as a boundary condition. Therefore, the plasma᎐air boundary velocity and the shock wave pressure in air are calculated self-consistently. We also assumed that the heat flux vanishes at the plasma᎐air boundary. We have checked explicitly, using complementary Hugoniot relations, that, in the early stage of the plasma expansion, the temperature of the shock wave can be high enough Žthat is, several tens of electron-volts. to dissociate and even ionize the air molecules, even for the modest laser fluences considered in this paper Ži.e. 10 Jrcm2 or smaller., a result which contradicts the weak shock wave assumption made here. The development of a better equilibrium shock wave model taking into account the phenomena of dissociation and ionization is currently in progress. Fortunately, the maximum plasma expansion predicted by the present model is quite reasonable, as discussed below.

3. Results 3.1. Laser ablation In this section, we present simulation results for the total amount of ablated matter, and for

the time characterizing the ablation, as a function of the laser fluence, and we compare them with the available experimental results. Although the process of mass removal considered here occurs on a much shorter time scale than the time scale of interest for the spectroscopy analysis used in LIPS Žby a factor smaller than approx. 10y3 ., its relevance for this paper is that it is an inherent feature of the phenomenon of plasma formation, and that it can be used for complementary checking of the validity of the code. The mechanisms of laser ablation shown by the code are given here in outline only because of their relative complexity. A more detailed discussion is presented elsewhere w27x. Initially, the laser pulse deposits energy at the solid surface over a distance of approximately the optical skin depth Žapprox. 10 nm for aluminum, for wavelengths of approx. 1 ␮m w28x.. This local heating drives a local rise in hydrodynamic pressure, which causes both the ejection of hot matter away from the target and the launching of a shock wave into the solid w9x. The code indicates that a clear separation between the ejected matter and the surviving target appears at a welldefined time Žhereafter called ablation time.. This result is the consequence of the fact that each homogeneous Lagrangian cell of the code can follow the Van der Waals type of behavior Žpresented by QEOS., including superheating and supercooling w29x. A close examination of the code results shows that the process of separation takes place through spinodal decomposition w30x, which implies that some matter elements enter the unstable zone Žnear the critical point of the phase diagram., where Ž⭸ pr⭸␳ .T ) 0. In the unstable zone, the least perturbation will cause a transition of the matter elements toward one of the two pressure extrema Ž⭸ pr⭸␳ .T s 0 boarding the unstable zone, i.e. toward either the vapor phase or the liquid phase. The end of the separation process is generally preceded by the formation of a series of large density oscillations, which can naturally be associated with ‘droplets’ and ‘bubbles’. The formation mechanism of these oscillations is of the same nature as the one described above for the separation between the remaining solid and the ablated

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

matter. The nearly solid-density ‘droplets’, prove to be permanent structures, which are ejected from the target with the rest of the ablated matter Žbeing in the vapor phase.. The existence of dense fragments in the ablated matter for ultrashort laser pulses has been suggested recently as an interpretation of experiments where Newton interference rings were observed nanoseconds after the laser impact w31x. ŽDespite this physically suggestive behavior, it should not be forgotten, however, that the model presented here lacks many features to provide an accurate description of an emulsion made of real droplets and vapor bubbles. In particular, the code cannot take into account the phenomena of evaporation and condensation that could change surface masses by transfer from one cell to another.. Fig. 1 shows the ablation depth and ablation time as a function of the incident laser fluence. Typically, the ablation depth is some tenths of micrometers while the ablation time is some tenths of nanoseconds. The reflection coefficient for the incident laser energy shows the transition from a metal-like behavior Žwhere the reflection coefficient increases as a function of the fluence. to a plasma-like behavior Žwhere the reflection coefficient decreases as a function of the fluence. observed experimentally w4᎐7x. The calculated ablation threshold is approximately 0.8 Jrcm2 , which is of the same order of magnitude as the values measured for several metals w32x. The ablation depth as a function of the fluence is also in the range of values measured for various metals w32,33x. It is likely, however, that the threshold ablation fluence is here somewhat overestimated and the ablation depth somewhat underestimated, since the equation of state used Ži.e. QEOS. probably does not describe exactly the thermodynamic properties of aluminum. In particular, the temperature at the critical point where all phases merge is likely overestimated by a factor of approximately two Žit should likely be 0.55 eV w34x instead of the 1.04-eV given by QEOS., so that more laser energy is required in this model to vaporize aluminum than is actually needed. While experiments on copper w33x for fluences

977

Fig. 1. Ablation depth and time when the ablation ends Žleft axis., and reflection coefficient Žright axis. as a function of the incident laser fluence. The laser parameters used here are: wavelength s 1 ␮m, pulse length s 500 fs, and normal angle of incidence.

in the range 1᎐10 Jrcm2 exhibit the linear dependence of the ablation depth as a function of the logarithm of the fluence, our results have significant curvature, predicting a rather more rapid increase. It is true, however, that, over a similar range of fluence as for that experiment, one could make a plausible straight-line fit in that range Žsee Fig. 1. and so claim agreement in this range. 3.2. Plasma expansion As shown in Fig. 1, the maximum duration of the ablation process, occurring for the largest fluence considered Ži.e. 100 Jrcm2 . is approximately 300 ps. We now present simulation results on much larger time scales, ranging between 0.1 and 20 ␮s. This range is essentially dictated by the needs for LIBS analysis. Indeed, it has been found experimentally that, for times earlier than approximately 0.1 ␮s, the signal-to-noise ratio is too small to provide any useful information about the density and the temperature of the plasma, since these quantities are inferred from the analysis of the spectral lines. For longer times Ži.e. for times beyond something like 10 ␮s, depending on

978

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

the species considered., the radiation emission decreases and eventually becomes too small to be useful. In order to simulate the plasma expansion over tens of microseconds for the least computer time, we restarted the simulation some time after the ablation time Žshown in Fig. 1., but only for the ablated matter, i.e. without the remaining solid target. All the quantities characterizing the ablated matter Ždensity, temperature, and velocity. were kept unchanged at the beginning of the simulation. The position of the interface between the plasma and the remaining bulk of solid was reset to z s 0. We assumed that, at this point, the velocity, the pressure and the heat flux all vanish. We checked that these simplification are of little consequence for the results for the expanding plasma presented in this section since, after the ablation process ends, there is little interaction between the remaining solid target and the expanding plasma. In the following, we present simulation results for a pulse length of 100 fs, a wavelength of 0.8 ␮m and a fluence of 10 Jrcm2 . Fig. 2 shows the position of the interface between the ablated matter and the air as a function of time for three different fluences. One can see that the interface position increases regularly as a function of time, following, initially, the approximate relation zA t 1r2 , and saturates at times that increase as the fluence increases. The saturation times can be estimated as 1, 3 and 7 ␮s for the fluences 2.5, 5 and 10 Jrcm2 , respectively. The maximum extension length reached by the plasma for a fluence of 10 Jrcm2 is approximately 2 mm, which is of the same order of magnitude as the measurements made for longer laser pulses with similar fluences w35x. The interface velocity gradually decreases as a function of time and finally vanishes at about the saturation times mentioned above, as a consequence of the resistance exerted on the plasma by the ambient air. For the fluence of 10 Jrcm2 , the plasma extension is larger than the measured laser spot diameter on the target of 0.6 mm. This implies that the use of the Cartesian geometry made here would likely fail at late times. A better approximation

Fig. 2. Position of the interface between the ablated matter and the air as a function of time for different fluences.

would consist of using the spherical geometry and assuming that the plasma expansion takes place within a given cone intersecting the laser heated spot on the target. However, in absence of a fully two-dimensional code, the angle ⍀ of that cone would have to be set rather arbitrarily. For that reason, we have chosen to present results only for the specific case where the angle ⍀ vanishes Ži.e. for the Cartesian geometry.. Figs. 3᎐5 show the evolution of the axial profiles of the atom and electron number density and of the temperature, respectively, at various times after the laser pulse, for a laser fluence of 10 Jrcm2 . These three quantities form the whole basis for the comparisons with the experiments. The axial profiles shown on these figures stop at positions corresponding to the extension of the plasma at the given time, as shown in Fig. 2. Note that all the axial profiles presented below ŽFigs. 3᎐7. do not start exactly at zs 0. The explanation for this is that the large density oscillations Ži.e. the ‘droplets’ and ‘bubbles’ making part of the ejected matter., mentioned in Section 3.1, have been erased because they would hide much of the more useful results. The profiles actually start right after the farthest ‘droplet’. ŽSince the ‘droplets’ are ejected from the remain-

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

979

Fig. 3. Atom number density profiles at different times after the laser pulse. The laser fluence is 10 Jrcm2 . The position z s 0 corresponds to the interface with the remaining target.

Fig. 5. Temperature profiles of the expanding plasma at various times. The laser fluence is 10 Jrcm2 .

ing solid target, the position of the farthest ‘droplet’ with respect to the position zs 0 increases as a function of time.. In any case, as discussed below, the matter to the left of this region does not contribute significantly to the plasma radiation emission.

One observes in Fig. 3 that the atom number density varies by a factor of approximately 10 or more over the axial profiles for all the times considered, even at 20 ␮s, where the plasma expansion has almost stopped Žsee Fig. 2..

Fig. 4. Electron density profiles of the expanding plasma at various times. The laser fluence is 10 Jrcm2 .

Fig. 6. Total pressure profiles of the expanding plasma at various times. The laser fluence is 10 Jrcm2 .

980

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

Fig. 7. Intensity of the 373.48-nm Fe I line in the aluminum plasma as a function of the axial position for various times. The laser fluence is 10 Jrcm2 .

The electron density shown in Fig. 4 is obtained by multiplying the average degree of ionization per atom ² Z : given by the code Žmore specifically, by the Thomas᎐Fermi ionization model, incorporated in QEOS. with the atom number density ŽFig. 3.. The electron number density profiles vary by a factor of 5 at most for the earlier time considered Ž100 ns. over the whole axial profiles, and by a factor of 2 only for the latest times considered. Since ² Z : s n ern a , one can see, using Figs. 3 and 4, that the concentration of ions is more important in the outer plasma. Fig. 5 shows that the plasma temperature globally decreases as a function of time. This was expected since, after the laser pulse, the plasma can only cool down as no other source of energy was supplied to the plasma. In our simulations, the plasma cooling is essentially due to the work of expansion of the plasma including that done on the ambient air Žon the right boundary. since we have not included in the code other loss mechanisms, such as radiation losses or diffusion Žthermal or particle. between the aluminum plasma and the colder ambient air. Within the range of time considered in Fig. 5, the maximum temperature in the profiles decreases only by approximately 4000 K. However, the simulation results

for earlier times Žnot shown here. indicate that, within the first 0.1 ␮s, the maximum temperature drops from approximately 10 5 K to 10 4 K. Using the results of Figs. 3᎐5, one can calculate the total pressure using the relation Žwhich is valid when the atom density is small enough., p s Ž na q ne . k BT

Ž6.

where n a and n e are the atom and electron number density, respectively. Fig. 6 shows that the pressure decreases and becomes more uniform as a function of time. At 20 ␮s, the plasma is practically in pressure equilibrium at the atmospheric pressure over the whole axial profile. This result is consistent with Fig. 2, since, for the fluence of 10 Jrcm2 , the plasma expansion has practically stopped at 20 ␮s. However, even at 20 ␮s, the plasma is still far from equilibrium with the ambient air since its maximum temperature is approximately 20 times the ambient air temperature Žsee Fig. 5. and the particle density Žatoms and electrons. at the plasma edges approximately 20 times smaller than the molecule density in ambient conditions Žsee Figs. 3 and 4.. ŽThe differences of shape in the pressure profiles observed for t - 5 ␮s are due to weak shock waves bouncing back and forth within the plasma during its

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

expansion. Such effects have already been noticed in w14᎐16x.. In order to compare our simulation results with our space-integrated and time-resolved experiments, we have also simulated the process of measurement by integrating the intensity emitted by each portion of the plasma and then extracted from this integrated signal an average temperature and an average density, using the same methods as the ones used for the experimental data. ŽObviously, this procedure is used here only because the plasma is inhomogeneous, otherwise it would be simpler to use directly the uniform temperature and electron density.. The excitation temperature has been determined from the experimental data using the Boltzmann method for a set of eight spectral lines of Fe Žan element that is initially present with a concentration of 0.7% in the aluminum target used.. The electron density has been determined using the usual Stark broadening formula w36x for the 281.62-nm line of Al II and the 285.21-nm line of Mg I ŽMg being, as Fe, an element that can be found with a concentration of 0.028% in the aluminum target used.,

˚. ⌬␭ 1r2 f 2Wn er10 16 Ž A

Ž7.

where ⌬␭ 1r2 is the full width at half maximum of the Stark broadened line, n e is the electron density in cmy3 and W is the electron impact ˚ We refer the reader to Le Droparameter in A. goff et al. w3x for a complete description of the experimental procedure. To calculate the intensity Žpower emitted per unit volume. for a transition from level j to level i, we used the formula: ␧ ji Ž T . s

hc n Ž T . g A eyE j r k B T 2␭ ji Q Ž T . j ji

Ž8.

Žwhere we have integrated over a solid angle of 2 ␲ .. In this expression, n is the number density of the emitting species Žneutral atoms or ions., ␭ ji is the wavelength of the radiation emitted, Q is the partition function, g j is the degeneracy factor and Ej is the excitation energy of level j.

981

As mentioned previously, we know from the experiments that for times greater than approximately 0.1 ␮s, some spectral lines of interest Žsuch as Al II lines Fe I and Mg I lines. emerge from the background. It has been checked carefully that these measured lines present no sign of auto-absorption, indicating that the plasma is transparent beyond that time. Fig. 7 shows the intensity profiles of the 373.48nm Fe I line, as obtained from Eq. Ž8., for different times. Unfortunately, the Thomas᎐Fermi ionization model used in the code provides only the average degree of ionization ² Z : Žfor a given atom density and temperature. and gives no information about the number density n appearing in Eq. Ž8.. The only exception occurs when ² Z : < 1, in which case one can see that the singly ionized atom density nq is almost equal to n e . However, by comparing the number density of atoms and of electrons in Figs. 3 and 4, respectively, one observes that ² Z : f 1 in the outer plasma, so that the relation nqs n e does not seem to be appropriate here, since the neutral atom density would be negligible. To overcome this difficulty, we have estimated n by means of the Saha equations, using the code atom number density ŽFig. 3. and plasma temperature ŽFig. 5.. However, as will be shown below, the Thomas᎐Fermi electron density can differ significantly from the Saha electron density since the physical basis of the two models are quite different. ŽWhile the Saha model can be considered more rigorous and more accurate than the Thomas᎐Fermi model at low density and low temperature, only the latter model can be used at nearly solid density.. For that reason, Fig. 7 is thus not fully consistent with the simulation results and must be considered as an indication only. The intensity profiles shown in Fig. 7 are very similar to the temperature profiles of Eq. Ž5. since the intensity ␧ ji wEq. Ž8.x is very sensitive to this variable while the atom density profiles ŽFig. 2. are too smooth to make significant changes. Since, as seen in Fig. 7, the most emissive zone is the outer region of the plasma Žthat corresponds to the earlier and hotter ablated matter., one infers that the measured plasma properties are actually the signature of this region of the plasma.

982

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

We have checked that the same conclusion holds for transition wavelengths different than the one considered in Fig. 7. One must remember, however, that for times as long as 20 ␮s, some mixing of the ablated matter with the colder ambient air could have been achieved Ža phenomenon sometimes called plasma collapse.. It is then likely that, in that a transition layer made of a mixture of ablated matter and of air molecules, the aluminum temperature in the outer plasma would be smaller than predicted here. It is worth noting that the most inner region of the plasma, i.e. the ‘bubbles-and-droplets’ layer Žsee Section 3.1. do not contribute to the emission for the transitions of interest here since the metallic ‘droplets’ are almost at the solid density, so that the value of the average degree of ionization ² Z : is high Ži.e. approx. 3 for aluminum.. For such values of ² Z :, the photon mean-freepath is very small since inverse Bremsstrahlung absorption is very significant. To calculate the code average density and temperature to be compared with the experiments, we integrated the intensity of a given spectral line ␧ ji Ž z . was given by Eq. Ž8.x of width ⌬␭ 1r2 Ž z . was given by Eq. Ž7.x over the axial positions z. A Lorentzian spectral line shape has been chosen for simplicity. The integrated intensity Žnow a surface power emission. has been used to calculate the average Boltzmann temperature, and the average electron density was obtained by using again Eq. Ž7.. The code average temperature and electron density are shown in Figs. 8 and 9 respectively, in comparison with the experimental results. Fig. 8 shows the average plasma temperature obtained from the code as well as the experimental results for the excitation temperature ŽTexc . and the ionization temperature ŽTioniz .. ŽThe ionization temperature has been determined by using the ratio between the Mg I line at 285.21 nm and the Mg II line at 279.55 nm.. By comparing the simulation results of Fig. 8 with Fig. 5, one can see that the average temperature is accurately given by the temperature at the outer plasma, which is as expected, since the inner part of the plasma contributes negligibly to the emission. One observes in Fig. 8 that the agreement

Fig. 8. Average plasma temperatures as a function of time for a laser fluence of 10 Jrcm2 . Texc and Tioniz refer, respectively, to the excitation and ionization temperature. The uncertainties on the experimental results are estimated as approximately 20%. w3x

between the simulation and the experiment here is very good since the calculated temperature lies between the two experimental temperatures. Fig. 9 shows that the agreement between the measured electron density and that predicted by the code Žin particular with the Thomas᎐Fermi ionization model. is not as good as for the temperature ŽFig. 8.. While the agreement can be judged as reasonable for short times Ž t - 0.3 ␮s., the experimental electron density decreases much faster than the simulation results and the discrepancy increases as a function of time. The reasons for this discrepancy are partly due to the use of the Thomas᎐Fermi ionization model, which, as mentioned before, is not very accurate at low density and low temperature. To illustrate this point, we have also shown in Fig. 9 the average electron density obtained by solving the Saha equations, using the code temperature and atom number density. The average electron density has then been obtained from the signal integrated over the axial positions, as explained above. One observes in Fig. 9 that the Saha electron density is now closer to the experimental results

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

Fig. 9. Average plasma density as a function of time for a laser fluence of 10 Jrcm2 . The uncertainties on the method are estimated as approximately 20%. w33x

than the Thomas᎐Fermi electron density although the measurements are still smaller. At approximately 4 ␮s, the Saha electron density is approximately five times the experimental density. The reasons for the discrepancies between the experimental and simulated electron densities are not fully clear yet. A discussion of the approximations made in the model and of their possible improvements will be presented in the next section. Finally, in order to check the consistency of the experimental results with the Saha equations, i.e. with the condition of thermodynamic equilibrium, we have calculated the plasma pressure, using the measured temperature and electron density. If the measurements were consistent with thermodynamic equilibrium, one would naturally expect the plasma pressure to be not too far from the atmospheric pressure a few microseconds after the laser pulse. Solving the Saha equations for the ion number densities nŽi ky1., where k y 1 is the degree of ionization Ž k - 4 proved to be ample enough., the atom density is simply n a s Ý ks 1 nŽi ky1. and the total pressure is obtained using Eq. Ž6.. For completeness, we have also shown the code average pressure using the code

983

average temperature ŽFig. 8. and density ŽFig. 9 ᎏ the Thomas᎐Fermi result.. One observes in Fig. 10 that the code average pressure tends towards atmospheric pressure as time increases, in agreement with Fig. 6. The average pressure inferred from the experimental electron density and excitation temperature is relatively close to the atmospheric pressure with reasonable precision. ŽThe excitation temperature was used because its calculation does not involve the electron density.. The fact that the pressure goes as low as 0.4 atm Žwhich, physically, does not make sense. can be due to experimental uncertainties. This seems to confirm that the condition of thermodynamic equilibrium is respected in the experiment with a reasonable precision, in agreement with one of the basic assumptions of our model. Actually, this conclusion was expected since: Ž1. as shown in Fig. 8, the measured excitation and ionization temperatures are very similar; and Ž2. the well-known necessary condition for thermodynamic equilibrium w37x: 3 n e G 1.6= 10 12 T 1r2 Ž ⌬ E . Ž cmy3 .

Ž9.

Fig. 10. Average pressure as a function of time. The ‘experimental’ pressure was obtained by solving the Saha equations, using the experimental excitation temperature and electron density. The laser fluence is 10 Jrcm2 .

984

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

Žwhere T is in Kelvin and ⌬ E in eV. is satisfied for all our measurements for ⌬ E- 4.5 eV, the value of 4.5 eV corresponding to the largest energy transition considered in our experiments w3x.

4. Summary and conclusion In this paper, we have presented simulation results for ultrashort laser pulse ablation and the subsequent plasma expansion in ambient air, and have compared them with the available experimental results. The one-dimensional fluid model used includes a self-consistent treatment of hydrodynamics, laser energy absorption, thermal conduction, and thermodynamics properties. Important approximations of the model are: Ži. the plasma is always in thermodynamic equilibrium; Žii. the radial plasma expansion is negligible; Žiii. the shock wave in air can be described using an Hugoniot shock wave model for cold air; and Živ. other energy loss mechanisms Žsuch as radiation emission and diffusion with the ambient air. are negligible. An important feature of the model discussed in this paper is that it involves no adjustable parameters. The results of the code for the ablation depth and the ablation threshold fluence appear to be in good agreement with measurements for metals. The code predicts the ejection of dense fragments shortly before the end of the ablation process Žwhich lasts typically hundreds of picoseconds., a result that also seems to be in agreement with experiments. The plasma expansion features have been calculated for times ranging between 0.1 and 20 ␮s and for a fluence of 10 Jrcm2 . The simulations show that for any of the times considered, the temperature and atom density can vary by a factor of approximately 10 within the axial plasma profile, while the electron density profiles are somewhat flatter. The plasma spectral emission profiles were found to be very similar to the temperature profiles and present a maximum in the outer plasma. From these emission profiles, an emission average temperature and electron

density were obtained for comparison with the experimental results. The code average temperature is in good agreement with experiments in most of the range of time considered Žexcept perhaps for times larger than 10 ␮s where the code temperature seems to decrease more slowly than in experiments.. However, the agreement is not as good for the electron density. The simulation results, using the Thomas᎐Fermi ionization model, are a factor of approximately six larger than the measured electron density at 100 ns, while this factor increases to 30 at 4 ␮s. We have shown, however, that these discrepancies can be partly explained by the fact that the Thomas᎐ Fermi model used in the code is not very accurate for the lower atom densities and lower temperatures. An estimate of the electron density made by means of the Saha equations Žusing the code atom density and temperature. reduces the discrepancy between the simulations and the experiments to a factor less than five. Having in view the complexity and the variety of the phenomena to simulate, and that the model contains no adjustable parameters, the performances of the model can be judged as rather satisfactory. However, in order to improve its accuracy, it will be necessary to revise and improve some of the approximations and simplifications made. In particular, when the plasma extension is larger than the laser spot size, it is likely that the radial expansion is important. As mentioned in Section 3.2, this effect can be easily taken into account by using the spherical geometry Žinstead of the Cartesian geometry. and assuming that the expansion takes place within a cone of angle ⍀ Ža parameter that needs to be set more or less arbitrarily .. The other possible improvements that could be brought to the model are not so straightforward, however. Diffusion between the ablated matter and the ambient air, radiation losses, and modeling of shock waves in air, taking into account dissociation and ionization, are important topics which are currently being considered. The forthcoming space-resolved measurements will be of great interest for checking more fully the code results.

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

Acknowledgements This work was financed partly by the National Science and Engineering Research Council of Canada ŽNSERC. under Strategic Grant STR 21898-98. References w1x D.A. Rusak, B.C. Castle, B.W. Smith, J.D. Winefordner, Fundamentals and applications of laser-induced breakdown spectroscopy, Crit. Rev. Anal. Chem. 27 Ž1997. 257᎐290. w2x V. Majidi, M.R. Joseph, Spectroscopic applications of laser-induced plasmas, Crit. Rev. Anal. Chem. 23 Ž1992. 143᎐162. w3x B. Le Drogoff, J. Margot, M. Chaker, M. Sabsabi, O. Barthelemy, T.W. Johnston, S. Laville, F. Vidal, Y. von ´ Kaenel, Temporal characterization of femtosecond laser pulses induced plasma for spectrochemical analysis of aluminum alloys, Spectrochim. Acta Part B 56 Ž2001. 987᎐1002. w4x A. Ng, P. Celliers, A. Forsman, R.M. More, Y.T. Lee, F. Perrot, M.W.C. Dharma-wardana, G.A. Rinker, Reflectivity of intense femtosecond laser pulses from a simple metal, Phys. Rev. Lett. 72 Ž1994. 3351᎐3354. w5x D.F. Price, R.M. More, R.S. Walling, G. Guethlein, R.L. Shepherd, R.E. Stewart, W.E. White, Absorption of ultrashort laser pulses by solid targets heated rapidly to temperatures 1᎐1000 eV, Phys. Rev. Lett. 75 Ž1995. 252᎐255. w6x K. Eidmann, J. Meyer-ter-Vehn, T. Schlegel, S. Huller, ¨ Hydrodynamic simulation of subpicosecond laser interaction with solid-density matter, Phys. Rev. E 62 Ž2000. 1202᎐1214. w7x M. Borghesi, A.J. Mackinnon, R. Gaillard, O. Willi, Absorption of subpicosecond UV laser pulses during interaction with solid targets, Phys. Rev. E 60 Ž1999. 7374᎐7381. w8x P. Celliers, A. Ng, Optical probing of hot expanded states produced by shock release, Phys. Rev. E 47 Ž1993. 3547᎐3565. w9x A. Ng, A. Forsman, P. Celliers, Heat front propagation in femtosecond-laser-heated solids, Phys. Rev. E 51 Ž1995. R5208᎐R5211. w10x A. Ng, P. Celliers, G. Xu, A. Forsman, Electron-ion equilibration in a strongly coupled plasma, Phys. Rev. E 52 Ž1995. 4299᎐4310. w11x R. Evans, A.D. Badger, F. Fallies, ` M. Mahdieh, T.A. Hall, P. Audebert, J.-P. Geindre, J.-C. Gauthier, A. Mysyrowicz, G. Grillon, A. Antonetti, Time-and-spaceresolved optical probing of femtosecond-laser-driven shock waves in aluminum, Phys. Rev. Lett. 77 Ž1996. 3359᎐3362.

985

w12x A.M. Komashko, M.D. Feit, A.M. Rubenchik, M.D. Perry, P.S. Banks, Simulation of material removal efficiency with ultrashort laser pulses, Appl. Phys., A 69 ŽSuppl.. Ž1999. S95᎐S98. w13x M.D. Feit, A.M. Rubenchik, B.-M. Kim, L.B. da Silva, M.D. Perry, Physical characterization of ultrashort laser pulse drilling of biological tissue, Appl. Surf. Sci. 127᎐129 Ž1998. 869᎐874. w14x J.N. Leboeuf, K.R. Chen, J.M. Donato, D.B. Geohegan, C.L. Liu, A.A. Puretzky, R.F. Wood, Modeling of dynamical processes in laser ablation, Appl. Surf. Sci. 96᎐98 Ž1996. 14᎐23. w15x J.N. Leboeuf, K.R. Chen, J.M. Donato, D.B. Geohegan, C.L. Liu, A.A. Puretzky, R.F. Wood, Modeling of plume dynamics in laser ablation processes for thin film deposition of materials, Phys. Plasmas 3 Ž1996. 2203᎐2209. w16x K.R. Chen, J.N. Leboeuf, R.F. Wood, D.B. Geohegan, J.M. Donato, C.L. Liu, A.A. Puretzky, Laser-solid interaction and dynamics of laser-ablated materials, Appl. Surf. Sci. 96᎐98 Ž1996. 45᎐49. w17x V.I. Mazhukin, I. Smurov, G. Flamant, Simulation of laser plasma dynamics: influence of ambient pressure and intensity of laser radiation, J. Comp. Phys. 112 Ž1994. 78᎐90. w18x V.I. Mazhukin, I. Smurov, G. Flamant, Simulation of laser plasma dynamics: influence of ambient pressure and intensity of laser radiation, Appl. Surf. Sci. 86 Ž1996. 303. w19x V.I. Mazhukin, S.G. Gorbachenko, I. Smurov, G. Flamant, A. Gleizes, 2D simulation of laser plasma dynamics and radiation transfer into the target, Proc. ECLAT-96 Ž1996. 603᎐612. w20x F. Garrelie, C. Champeaux, A. Catherinot, Study by a Monte Carlo simulation of a background gas on the expansion dynamics of a laser-induced plasma plume, Appl. Phys. A 69 Ž1999. 45᎐50. w21x W. Rozmus, V.T. Tikhonchuk, A model of ultrashort laser pulse absorption in solid targets, Phys. Plasmas 3 Ž1996. 360᎐367. w22x J.P. Christiansen, D.E.T.F. Ashby, K.V. Roberts, MEDUSA a one-dimensional laser fusion code, Comput. Phys. Commun. 7 Ž1974. 271᎐287. w23x Y.T. Lee, R.M. More, An electron conductivity model for dense plasmas, Phys. Fluids 27 Ž1984. 1273᎐1286. w24x L. Ginzburg, The Propagation of Electromagnetic Waves in Plasmas, Wiley, New York, 1970. w25x R.M. More, K.H. Warren, D.A. Young, G.B. Zimmerman, A new quotidian equation of state ŽQEOS. for hot dense matter, Phys. Fluids 31 Ž1988. 3059᎐3078. w26x Ya.B. Zel’dovich, Yu.P. Raizer, Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena, section 15, Academic Press, New York, 1966, p. 50 w27x F. Vidal, T.W. Johnston, S. Laville, O. Barthelemy, M. ´ Chaker, B. Le Drogoff, J. Margot, M. Sabsabi, Criticalpoint phase separation in laser ablation of conductors, Phys. Rev. Lett. 86 Ž2001. 2573᎐2576.

986

F. Vidal et al. r Spectrochimica Acta Part B: Atomic Spectroscopy 56 (2001) 973᎐986

w28x M. von Allmen, A. Blatter, Laser-beam interaction with materials. Physical principles and applications. Springer Ser. Materials Sci., 2nd ed., Springer-Verlag, Berlin, Heidelberg, 1995. w29x L.D. Landau, E.M. Lifshits, Course of Theoretical Physics V: Statistical Physics, chapter 8, Pergamon Press, Oxford, 1978. w30x V.P. Skripov, A.V. Skripov, Spinodal decomposition Žphase transitions via unstable states ., Sov. Phys. Usp. 22 Ž1979. 389᎐410. w31x K. Sokolowski-Tinten, J. Bialkowski, A. Cavalleri, D. von der Linde, A. Oparin, J. Meyer-ter-Vehn, S.I. Anisimov, Transient states of matter during short pulse laser ablation, Phys. Rev. Lett. 81 Ž1998. 224᎐227. w32x S. Preuss, A. Demchuk, M. Stuke, Sub-picosecond UV laser ablation of metals, Appl. Phys. A 61 Ž1995. 33᎐37. w33x S. Nolte, C. Momma, H. Jacobs, A. Tunnermann, B.N. ¨

w34x

w35x

w36x w37x

Chichkov, B. Wellegehausen, H. Welling, Ablation of metals by ultrashort laser pulses, J. Opt. Soc. Am. B 14 Ž1997. 2716᎐2722. A.V. Bushman, I.V. Lomonosov, V.E. Fortov, Equations of State of Metals at High Energy Densities, Institute of Chemical Physics, Russian Academy of Sciences, Chernogolovka, 1992 in Russian. K.J. Grant, G.L. Paul, Electron temperature and density profiles of excimer laser-induced plasmas, Appl. Spectroscopy 44 Ž1990. 1349᎐1354. H.R. Griem, Plasma Spectroscopy, McGraw-Hill, New York, 1964. R.W.P. McWhirter, Spectral Intensities in: R.H. Huddlestone, S.L. Leonard ŽEds.., Plasma Diagnostic Techniques, Chapter 5, Academic Press, New York, 1965, p. 206