Eighth Mississippi State - UAB Conference on Differential Equations and Computational Simulations. Electronic Journal of Differential Equations, Conference 19 (2010), pp. 1–14. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu ENTHALPY MODEL FOR HEATING, MELTING, AND VAPORIZATION IN LASER ABLATION VASILIOS ALEXIADES, DAVID AUTRIQUE ABSTRACT. Laser ablation is used in a growing number of applications in various areas including medicine, archaeology, chemistry, environmental and materials sciences. In this work the heat transfer and phase change phenomena during nanosecond laser ablation of a copper (Cu) target in a helium (He) background gas at atmospheric pressure are presented. An enthalpy model is outlined, which accounts for heating, melting, and vaporization of the target. As far as we know, this is the first model that connects the thermodynamics and underlying kinetics of this challenging phase change problem in a self-consistent way. 2000 Mathematics Subject Classification. 92C45, 35K60, 65M99. Key words and phrases. Laser ablation; phase change; enthalpy scheme; thermochemistry; finite volume scheme. © 2010 Texas State University - San Marcos. Published September 25, 2010. 1. INTRODUCTION Laser ablation (LA) is the process of removing material from a target using a laser beam. It is used in a growing number of applications, such as pulsed laser deposition, nanoparticle manufacturing, micromachining, surgery, as well as chemical analysis [30, 32, 11]. Direct solid microanalysis using LA in combination with inductively coupled plasma mass spectrometry (LA-ICP-MS) is nowadays one of the most frequently used techniques for a fast and powerful multi-element determination of solid samples at the trace and ultra-trace concentration levels of a wide variety of sample types [31, 8, 4, 14]. The growing interest in LA as a sample introduction technique stems from the ability to sample diverse materials, ranging from conducting to non-conducting, inorganic and organic compounds as solids or powders. Besides bulk analysis, the focussing characteristics of lasers permit sampling in small areas, so that localized microanalysis and spatially resolved studies become feasible. The processes occuring during LA-ICP-MS are depicted schematically in Fig. 1. FIGURE 1. Schematic of the Laser Ablation - Inductively Coupled Plasma Mass Spectrometry (LA-ICP-MS) process. (Courtesy of C.C. Garcia, ISAS, Dortmund) Within an ablation cell, an amount of matter is extracted from the sample by focusing a laser beam on its surface. The sample heats up, melts, vaporizes, and expands against a carrier gas at atmospheric pressure. The expanding plume partially ionizes under the laser radiation and a plasma is generated above the target. Complex physical processes occuring at the target surface and in the expanding plume result in the formation of small particles with sizes ranging from several nanometers to a few micrometers. The aerosol is subsequently transported by the carrier gas through a tubing system to an ICP-torch where it is ionized. Finally a small fraction of the ions is sent to the MS-instrumentation where they are analyzed. The laser ablation process is typically investigated by means of kinetic [28, 12, 18, 13, 24], hydrodynamic [7, 9, 15, 5, 26], or hybrid models [19, 20]. Because kinetic and hybrid models require quite long calculation times, one tries to apply hydrodynamic models if at all possible. Since the laser-target heating during LA-ICP-MS results in a cascade of complex physical processes (see §2), it is of utmost importance to include in the existing hydrodynamical models a consistent treatment of target heating. Here the kinetic relations that interconnect the physical processes in and above the target determine the pathway the metal follows on the phase diagram during laser heating. Accordingly, one has to account for them in order to arrive at a physically consistent treatment of laser-material interaction. In §2 we briefly describe some issues that arise during the modeling of nanosecond laser ablation. In §3 we describe a seamless enthalpy-based model treating the heating, melting and vaporization of the target. Numerical simulations for laser ablation of copper (Cu) in a background gas (He) in 1-D, are presented in §4, and final conclusions in §5. 2. NANOSECOND LASER ABLATION Several complex, tightly coupled physical processes occur in and just above the target. The target heats up, melts and vaporizes. Above the target, the evaporated particles achieve translational equilibrium within a few mean free paths by means of collisions. This transition layer is known as the Knudsen layer (KL). The vapor near the target is a partially ionized gas that is separated from the surface by an electrical sheath across which there is a potential drop that accelerates the charged particles. The Knudsen layer provides the connection between the ‘inner world’ (target) and the ‘outer world’ (plasma plume). Beyond this layer, the plasma absorbs laser energy before it can reach the target, giving rise to plasma shielding. The absorption of laser energy results in very high plume temperatures, velocities, species densities, pressures, and so on. At later times, small particulates are formed during the cooling stage of the plume. Here homogenous nucleation and recondensation result in nanosized particles. Moreover, recoil pressures acting on the melt can cause melt motion and melt expulsion forming larger particles. The melt pool dynamics leads to formation of jets at the edge of the formed crater, that break up in micro-droplets [32, 6]. An effective hydrodynamic model of the entire ablation process should posess the following highly desirable features: • seamless treatment of heating / melting / vaporization, • temperature-dependent thermophysical and optical properties, • treatment of vapor as a non-ideal gas by means of an equation of state(EOS) up to and above the critical temperature Tcrit , • treatment of the plasma-sheet at the target surface, • coupling target and plume by means of the Knudsen layer, • treatment of plume expansion and recondensation, in vacuum or in a background gas • dealing with plasma formation (ions and electrons) mechanisms, • treating absorption of laser energy in the plume and shielding of the target, • capturing (very) strong shocks via high resolution numerical schemes, • be applicable in 1D, 2D axisymmetric and 3D spatial dimensions. Additional challenges in modeling laser ablation arise from: • extreme space and time scales, • extreme gradients: temperature may rise to thousands of degrees locally, • extreme variation in thermophysical properties, • the need for extensive thermophysical data: T-dependent density, heat capacity, thermal conductivity, phase diagram for solid, liquid, vapor over the range 300 K to Tcrit ( ∼ 8000 K ), • the need for T-dependent and wavelength dependent optical data. 3. ENTHALPY MODEL FOR TARGET HEATING - MELTING - VAPORIZATION In this section we develop an enthalpy formulation for target heating - melting - vaporization which addresses the issues listed above. Heat conduction and phase changes in the target are naturally described by the energy conservation law in terms of the (volumetric) enthalpy H. ∂tH + ∇ · Q~ = S . (3.1) The laser source term S, at time t and position z, is given by S(t, z) = I(t)(1 − R) αe−αz , (3.2) with R and α the reflection and absorption coefficients of the target, respectively, and intensity I(t) = βIo e^{-(t-tmax)^2/2σ^2} . (3.3) Here β is the shielding coefficient at the target that depends on the plume temperature and species densities (neutrals, ions and electrons) above the target. The calculations presented later in §4 are performed for a Gaussian-shaped laser pulse of full width at half maximum (FWHM) 10 ns, and peak laser intensity of the order of 10^13 W/m^2, at a wavelength of 266 nm. In Fig. 3 the laser intensities before (β = 1) and after shielding at the surface are shown. FIGURE 2. Time course I(t) of laser pulse intensity with FWHM=10 ns, peak 7 × 10^12 W/m^2. There are three time scales that characterize laser heating: τe = Ce/γ, the electron cooling time, τi = Ci/γ, the lattice heating time, and τl, the duration of laser pulse. Here Ce and Ci are the heat capacities of the lattice and the electron subsystems, respectively, and γ a parameter characterizing the electron-lattice coupling in the target. Since τe is of the order of 1 ps and τl is of the order of ns, the electron absorbed energy has enough time to be transferred to the lattice [10]. Consequently, the electrons and the lattice reach thermal equilibrium and heat conduction in the target can be described by means of Fourier’s law provided that the heat flux is constant over several mean free paths of the electrons that carry the heat [17]. Since the melt thickness is very small, convective effects do not arise in the melt. Thus the heat flux in the solid and liquid is due to heat conduction only, Q~ = −k∇T, with k = thermal conductivity, T = temperature. On the contrary, the presence of strong advection in the plume domain requires solution of the continuity and momentum equations as well. Note that, depending on the local pressure P, vaporization can occur at any temperature above absolute zero, even below the melt point Tm as well as above the normal boiling point Tb. Therefore the vaporization temperature TV (P) is unknown and must be determined in the course of the solution. Finally note that the surface recession velocity, which determines the crater depth in the target, is also an unknown, since the surface is an unknown moving boundary. The energy conservation law (3) is valid mathematically, in the weak sense, throughout the system irrespective of phase (with coefficients appropriate to phase: k^S, k^L, k^V, etc.). It directly updates the enthalpy H (per unit volume) to new times, from which we need to find the new phase, temperature, and phase fractions. This requires Equations of State (EOS) H = H(T, P, phase) consistent with the thermochemistry of the material, which we describe below. 3.1. Phases from volumetric enthalpy. The development is based on the basic Gibbs relation for the exact differential of enthalpy [1] in any phase i = S, L, V : dH^i = C_p^i dT + [−α^i T + 1]dP i = S, L, V , (3.4) where C_p^i and α^i are the heat capacity (per unit volume) and thermal expansion coefficient of phase i, respectively. The enthalpy within each phase (Solid, Liquid, Vapor) can be expressed in terms of the pair (T, P) by integration along constant T and constant P paths on the (T, P)-phase diagram, of the Gibbs relation (3.1), relative to a consistent reference state. In what follows, upper case H denotes volumetric enthalpy (per unit volume), lower case h denotes specific enthalpy (per unit mass), and values with subscript "o" denote interfacial (switch) values. Clearly, H = ρh in each phase. We are dealing with three phases: solid (S), liquid (L), vapor (V). The amount of phase i = S, L, V present can be quantified by its volume phase fraction, Λ^i, or by its mass phase fraction, λ^i. The fractions are related by: ρ^i Λ^i = ρ λ^i (i = S, L, V), where ρ = total density (of mixture) = Σ_i ρ^i Λ^i. Note that fractions must add up to 1: Λ^S + Λ^L + Λ^V = 1 and λ^S + λ^L + λ^V = 1. For simplicity, we assume the density of solid, ρ^S, and of liquid, ρ^L, remain constant (but ρ^S ≠ ρ^L). As reference state we choose solid at temperature Tm and pressure Pref = 1 atm. The value href (and thus also Href = ρ^S href) of the enthalpy at the reference state can be arbitrary; it will be chosen later (so as to make ∆Hvap > 0). By integration of (3.1) within each phase, and changing phase at (Tm, Pref) (from solid to liquid) and at (TV, Pref) (from liquid to vapor) we can characterize the phases in terms of the value of H alone (which is being updated by the energy conservation law). This is done in terms of certain switch values H_o^S, H_o^L, H_o^{LV}(TV), and H_o^V(TV), defined below, with TV any desired vaporization temperature (Tm < TV ≤ Tb). • Solid phase (T ≤ Tm): H < H_o^S := Href ≡ ρ^S h_o^S, (set h_o^S := href). Then Λ^S = 1, Λ^L = 0, Λ^V = 0, ρ = ρ^S, H(T) = H_o^S + ∫_{Tm}^T C_p^S(τ) dτ. Thus, given H, we can find the temperature T (< Tm) by solving the equation: ∫_{Tm}^T C_p^S(τ) dτ = H − H_o^S for T. • Solid-Liquid transition (at T = Tm): H_o^S < H < H_o^L := ρ^L h_o^L, with h_o^L = h_o^S + ∆h^{fus}. Set ∆H^{fus} := H_o^L − H_o^S ≡ (ρ^L − ρ^S)h_o^S + ρ^L ∆h^{fus}, which can be made > 0 by choosing href to be ≤ 0. Then T = Tm, Λ^L = (H − H_o^L)/∆H^{fus}, Λ^S = 1 − Λ^L, Λ^V = 0, ρ = Λ^S ρ^S + Λ^L ρ^L, H = H_o^L + Λ^L ∆H^{fus}, h = H/ρ. • Liquid phase (Tm < T < TV): H_o^L < H < H_o^{LV}(TV) := H_o^L + ∫_{Tm}^{TV} C_p^L(τ) dτ. Then Λ^L = 1, Λ^S = 0, Λ^V = 0, ρ = ρ^L, H(T) = H_o^L + ∫_{Tm}^T C_p^L(τ) dτ. Thus, given H, we can find the temperature T (Tm < T < TV) by solving for T the equation: ∫_{Tm}^T C_p^L(τ) dτ = H − H_o^L. • Liquid-Vapor transition (at T = TV): H_o^{LV}(TV) < H < H_o^V(TV, Pref) := ρ^V(TV, Pref) h_o^V, with h_o^V(TV, Pref) := h_o^{LV}(TV) + ∆h^{vap}(TV, Pref), where ρ^V(TV, Pref) = vapor density at (TV, Pref), and ∆h^{vap}(TV, Pref) = heat of vaporization at (TV, Pref). Set ∆H^{vap} := H_o^V − H_o^{LV} ≡ (ρ^V − ρ^L)h_o^{LV}(TV) + ρ^V ∆h^{vap}(TV, Pref), which can be made > 0 for any Tm < T < Tcrit by choosing href sufficiently negative, namely href < −∆h^{fus} − c_p^L [Tcrit − Tm], thanks to the fact that ∆h^{vap}(Tcrit) = 0. Then Λ^V = (H − H_o^{LV})/∆H^{vap}, Λ^L = 1 − Λ^V, Λ^S = 0, ρ = Λ^L ρ^L + Λ^V ρ^V, H = ρh = H_o^{LV}(TV) + Λ^V ∆H^{vap}, h = H/ρ. • Vapor phase (T > TV): H > H_o^V(TV, Pref). Then Λ^V = 1, Λ^L = 0, Λ^S = 0, H(T, P) = H_o^{LV}(TV) + ∫_{TV}^T C_p^V(τ, Pref) dτ + ∫_{Pref}^P [−T α^V(T, p) + 1] dp, with α^V the thermal expansion coefficient of vapor. Since the heat capacity is strictly positive and ∆H^{fus} and ∆H^{vap} are also positive (by our choice of href), the dependence of H on T is monotonic, as shown schematically in Fig. 3.1. Thus T can be found from H. FIGURE 3. Schematic dependence of temperature T and phase fractions on volumetric enthalpy H. The PDE (3) is discretized by a finite volume method, and explicitly in time. Briefly, the enthalpy algorithm proceeds as follows: Knowing the values of T and phase fractions at time step tn, we evaluate the thermophysical property values, compute the flux Q~ , and update the enthalpy H from the discretized PDE to time tn+1. From H we find the phase indicators H_o^S, H_o^L, H_o^{LV}, H_o^V, which determine the temperature and phase fractions as described above. 3.2. Knudsen Layer and Vaporization temperature. Whereas melting occurs at a known melt temperature Tm (= 1358 K for copper), vaporization can occur at any temperature TV, up to the critical temperature Tcrit (∼ 8000 K for copper), and must be determined at each time step. TV depends on the pressure above the melt, which is coupled to the Knudsen layer. Across the Knudsen Layer, treated as a gas dynamic discontinuity (shock), conservation of mass, momentum and energy impose certain jump conditions on temperature, pressure, density, and velocity. There is extensive literature on modeling the Knudsen layer, at various scales and levels of detail, e.g. [2, 34, 22, 23, 27, 29, 12, 16, 15]. A useful treatment is given in Gusarov-Smurov [16]. Contrary to most earlier works, their treatment accounts for evaporation into a background gas. Moreover, it considers condensation back to the target, occuring when the outer pressure exceeds the surface pressure. The KL relations couple the target and plume, and the pressure ratio across the KL governs what happens at the melt surface. Let TV, PV be the values at the melt surface (vaporization front), and TKL, PKL denote the values at the outer side of the Knudsen layer (away from the target, towards the plume). If the pressure ratio Pratio := PKL/PV < 1 then subsonic evaporation occurs. If Pratio > 1 then condensation occurs (vaporized material re-condenses onto the melt). Note that in case of ablation into vacuum one can assume sonic evaporation, in which case the pressure and temperature ratios over the KL are constant, resulting in the decoupling of the target and plume processes. Every quantity experiencing a jump across the Knudsen layer is ultimately a function of the vaporization temperature TV, which changes in time and needs to be evaluated at each timestep. We have developed a new approach, based on enthalpy, that includes the underlying evaporation kinetics. The KL expressions from [16] are applied in an iterative procedure that determines TV in a vaporizing control volume. An outline of the approach (too involved to present here in detail): is as follows: The updated enthalpy H detects if a control volume is vaporizing (if H_o^{LV}(TV, Pref) < H < H_o^V(TV, Pref)) at the current TV. From the temperature TKL, pressure PKL, and gas density ρ^{KL} on the plume side, we find the recession velocity vV of the vaporization front (by the Hertz-Knudsen relation and conservation of mass). From the recession velocity vV we update the vapor fraction of the (control) volume Vol during a time-step ∆t via Λ^V = Λ_{old}^V + vV Vol/∆t. Then we estimate the new enthalpy of the volume as Hnew(TV) = H_o^{LV}(TV) + Λ^V ∆H^{vap}(TV), which should equal H. This constitutes an equation Hnew(TV) = H for the unknown TV, which we solve iteratively (by a bisection type root finder, such as Brent’s algorithm). Beyond the Knudsen layer, the plume of vapor expands rapidly and ionizes into plasma. The flow is modeled by the Euler equations of gas dynamics (plus species conservation for electrons, ions, neutrals, He gas), which must be solved simultaneously with the target since they are coupled at the Knudsen Layer. We employ finite volume discretization, explicit time-stepping, and a central 2nd order high resolution advection scheme. To save computational time, we use an adaptive grid in the plume which expands ahead of the plume. Our numerical scheme allows temperature and pressure dependence of all thermophysical properties, and vapor properties from equations of state (such as [25] for copper). 4. NUMERICAL SIMULATIONS We present 1-D simulations of laser ablation of Cu. The target is initially set at room temperature (300 K). Above the target we assume a stationary background gas (He), initially set at atmospheric pressure (1 atm) and room temperature (300 K). The entire sytem (target+background gas) is irradiated by a laser at the conditions mentioned in §3. The target, of thickness 12 µm, is discretized along the axial direction (direction of the laser beam), into a non-uniform grid of 2000 control volumes, denser near the surface. In the plume region we use (much) bigger cells, of length ∆z = 5 × 10^{−7} m = 500 nm, the number of which is adaptively increased to extend the grid beyond the expanding plume. The time-step is taken to be ∆t = 1 × 10^{−13} s, respecting the CFL condition in the entire computational domain. The total simulation time is 30 ns, which corresponds to a single pulse, cf. Fig. 3. Thermodynamic properties of Cu are obtained from the ITTEOS [25] EOS data (via interpolation). As a result of the deposited laser energy the target heats up very fast. The time evolution of the melting and vaporization fronts is shown in Fig. 4. At 7.5 ns melting starts and vaporization is assumed to be minor. At 10 ns the first cell exceeds the normal boiling point (Tb = 2835K) and vaporization becomes sigificant. Thus the vaporization (surface) temperature TV and pressure PV, as well as the pressure at the outer side of the KL, PKL, increase quickly as seen in Fig. 4–4. Since PKL < PV the evaporation is subsonic (see §3.2). FIGURE 4. Melt and vaporization depths vs time. When the surface reaches temperatures near the critical state, the melt becomes transparent. This bleaching of liquid properties occurs at temperatures about 0.9 Tcrit [3], [21], [33]. Here a strong reduction of the electron density in the liquid results in a drastic change of the optical properties. As the incident laser energy penetrates through this transparent layer to the underlying material, the transparency front starts to propagate into the interior liquid. The transparent surface cell now evaporates at a constant vaporization temperature TV ≈ 0.9 Tcrit as can be seen in Fig. 4. FIGURE 5. Vaporization temperature vs time. In this hot dense plume, an amount of electrons is created due to photoinization of the excited atoms. These electrons absorb additional energy in, so called, inverse bremssthralung processes [35]. Due to the aborption of laser energy, additional electrons and ions are formed resulting in an avalanche effect that keeps on increasing the plume temperatures to a few 10000 K (here 30000 K): a plasma is created. At a certain instant, here at ≈ 16 ns, the plasma is so dense that it completely shields the target. As a result, the vaporization temperature and pressure decrease, and when PV < PKL we arrive at the condensation regime (Fig. 4). During this stage, vapor from the plume recondenses on the target, the recession velocity becomes negative (Fig. 4), and the crater depth decreases (4). After a while, at ≈ 20 ns, the plume density above the target decreases due to advection, the laser light reaches again the target surface raising the surface temperature and pressure, and we return to the subsonic evaporation regime. Accordingly, the laser intensity profile at the target takes a bimodal shape (Fig. 3). As the applied laser pulse diminishes toward the end, the amount of laser energy deposited in the plume decreases, the surface temperature and pressure decrease as well, and material recondenses on the surface. In the final stage evaporation and condensation compete, resulting in a recession rate of ≈ 0 and constant evaporation depth of ≈ 150 nm, as seen in Fig. 4 – 4. FIGURE 6. Vaporization (saturation) pressure PV vs time, and pressure at the outer side of the Knudsen layer, PKL, for comparison. FIGURE 7. Vaporization (recession) velocity vs time. FIGURE 8. Vaporization depth vs time. Maximum depth reaches 200 nm shortly after 15 ns (when laser intensity peaks). 5. CONCLUSION A seamless enthalpy formulation of heat transfer and phase changes during nanosecond laser ablation of a target is presented. The formulation is based on Equations of State of the form H = H(T, P, phase) consistent with the thermochemistry of the material. It allows full temperature (and/or pressure) dependence of thermophysical and of optical properties, and can use available EOS data, up to critical temperature. Thus, it can realistically describe all the phases of real materials. The enthalpy model tracks both the solid/liquid and liquid/vapor fronts and accounts for the kinetic relations at the liquid/vapor front (Knudsen layer) using an iterative procedure. It couples naturally to plasma dynamics models in the plume and the coupled model posseses all the desirable features listed in §2. Thus, the entire laser ablation process can be simulated, from heating of the target, to plasma formation and laser absorption, to the shielding effect of plasma on the target and recondensation. As far as we know, this is the first model that connects the thermodynamics and underlying kinetics of this challenging phase change problem in a self-consistent manner. Being dimension-independent, it is applicable in 1D, 2D, and 3D modeling of laser ablation. As proof of principle, we presented a one-dimensional numerical implementation for laser ablation of a copper target in helium background gas at atmospheric pressure, coupled to a gas dynamics model in the plume (where electons, ions, and He gas are tracked). As we saw in §4, the expected physical behavior of the system is captured very well. Acknowledgements. We thank Drs. PR Levashov and KV Khishchenko for providing EOS data for Cu. We also thank C.C. Garcia of ISAS, Dortmund, for the drawing of the LA-ICP-MS-setup (Fig. 1). Support of the second author by the Flemish Fund for Scientific Research (FWO Vlaanderen) is gratefully acknowledged. REFERENCES [1] Vasilios Alexiades and Alan D. Solomon, Mathematical Modeling of Melting and Freezing Processes, Hemisphere Publ.Co., 1993. [2] S. I. Anismov, ”Vaporization of metal absorbing laser radiation”, Sov. Phys., JETP, 27: 182-183, 1968. [3] V. A. Batanov, F. V. Bunkin, A. M. Prokhorov, and V. B. Fedorov; “Evaporation of metals caused by intense optical radiation”, Sov. Phys., JETP, 36: 311-322, 1973. [4] J. S. Becker, “Applications of inductively coupled plasma mass spectrometry and laser ablation inductively coupled plasma mass spectrometry in materials science”, Spectrochim. Acta, Part B, 57: 1805–1820, 2002. [5] D. Bleiner, “Mathematical modelling of laser-induced particulate formation in direct solid microanalysis”, Spectrochim. Acta Part B, 60(1): 49-64, 2005. [6] D. Bleiner and A. Bogaerts, “Computer simulations of laser ablation sample introduction for plasma-source elemental microanalysis”, J. Anal. At. Spectrom., 21: 1161–1174, 2006. [7] Annemie Bogaerts, Zhaoyang Chen, Renaat Gijbels, Akos Vertes, “Laser ablation for analytical sampling: what can we learn from modeling?”, Spectrochimica Acta B, 58: 1867-1893, 2003. [8] J. A. C. Broekaert, Analytical Atomic Spectrometry with flames and plasmas, Wiley VCH, 2002. [9] Zhaoyang Chen, Annemie Bogaerts; “Laser ablation of Cu and plume expansion into 1 atm ambient gas”, J. Appl. Phys, 97: 063305, 2005. [10] B. N. Chichkov, C. Momma, S. Nolte, F. von Alvensleben, A. Tünnermann, “Femtosecond, picosecond and nanosecond laser ablation of solids”, Appl. Phys. A, 63: 109–115, 1996. [11] D. B. Chrisey and G. K. Hubler, Pulsed Laser Deposition of Thin Films, Wiley, 1994. [12] O. Ellegaard, J. Schou, H. M. Urbassek, “Monte Carlo description of gas flow from laser-evaporated silver”, Applied Physics A, 69: S577S581, 1999. [13] F. Garrelie, J. Aubreton, and A. Catherinot, “Monte Carlo simulation of the laser-induced plasma plume expansion under vacuum: Comparison with experiments”, J. Appl. Phys., 83, 5075, 1998. [14] D. Günther, “Laser ablation-inductively coupled plasma mass spectrometry trends”, Anal. Bioanal. Chem., 372: 31–32, 2002. [15] A. V. Gusarov, A. G. Gnedovets, and I. Smurov, “Two-dimensional gas-dynamic model of laser ablation in an ambient gas,” Applied Surface Science, 154-155: 66-72, 2000. [16] A. V. Gusarov, I. Smurov, “Thermal model of nanosecond pulsed laser ablation: Analysis of energy and mass transfer”, J. Applied Physics, 97: 014307, 2005. [17] R. E. Harrington, “Application of the Theory of Heat Conduction to the Absorption of Blackbody Radiation”, J. Applied Physics, 38: 3266, 1967. [18] T. E. Itina, V. N. Tokarev, W. Marine, and M. Autric, “Monte Carlo simulation study of the effects of nonequilibrium chemical reactions during pulsed laser desorption”, J. Chem. Phys., 106, 8905, 1997. [19] T. E. Itina, J. Hermann, Ph. Delaporte, and M. Sentis, “Laser-generated plasma plume expansion: Combined continuous-microscopic modeling”, Phys. Rev. E, 66, 066406, 2002. [20] T. E. Itina, J. Hermann, Ph. Delaporte, and M. Sentis, “Combined continuous–microscopic modeling of laser plume expansion”, Appl. Surf. Sci., 27: 208–209, 2003. [21] R. V. Karapetyan and A. A. Samokhin, “Influence of an increase in the transparency on the intense evaporation of metals by optical radiation”, Sov. J. Quantum Electron., 4: 1141–1142, 1975. [22] C. J. Knight, “Theoretical modeling of rapid surface vaporization with back pressure”, AIAA Journal, 17(5): 519–523, 1979. [23] C. J. Knight, “Transient Vaporization from a Surface into Vacuum”, AIAA Journal, 20(7): 950-954, 1981. [24] J. C. S. Kools, “Monte Carlo simulations of the transport of laser-ablated atoms in a diluted gas”, J. Appl. Phys., 74, 6401, 1993. [25] P. R. Levashov and K. V. Khishchenko, ITTEOS 5.8 software for calculation of EOS for metals, 2007. [26] H. Lindner, D. Autrique, C. C. Garcia, K. Niemax and A. Bogaerts, “Optimized transport setup for high repetition rate pulse-separated analysis in laser ablation-inductively coupled plasma mass spectrometry”, Anal. Chem., 81: 4241–4248, 2009. [27] A. Miotello, A. Peterlongo, R. Kelly, “Laser-pulse sputtering of aluminium: gas-dynamic effects with recondensation and reflection conditions at the Knudsen layer,” Nuclear Instruments and Methods in Physics Research B, 101: 148-155, 1995. [28] I. NoorBatcha, R. R. Lucchese, and Y. Zeiri, “Monte Carlo simulations of gas-phase collisions in rapid desorption of molecules from surfaces”, J. Chem. Phys., 86, 5816, 1987. [29] E. Ohmura, I. Fukumoto, I. Miyamoto, “Molecular dynamics simulation of laser ablation of metal and silicon”, Int. J. Jpn. Soc. Prec. Eng., 32: 248-253, 1998. Int. J. Jpn. Soc. Prec. Eng. 32 [30] J. F. Ready, Effects of High Power Laser Radiation, Academic Press, 1971. [31] A. Vertes, R. Gijbels, F. Adams, Laser Ionization Mass Analysis, Wiley VCH, 1993. [32] M. von Allmen, Laser-Beam Interactions with Materials, Springer, 1987. [33] J. H. Yoo, S. H. Jeong, X. L. Mao, R. Greif and R. E. Russo, “Evidence for phase-explosion and generation of large particles during high power nanosecond laser ablation of silicon”, Appl. Phys. Lett., 76: 783-785, 2000. [34] T. Ytrehus, “Theory and experiments on gas kinetics in evaporation”, pp.1197-1212 in Rarefied Gas Dynamics, edited by Potter, AIAA, 1977. [35] B. Zelldovich and Yu P. Raizer, Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena, Academic Press, 1966. Vasilios Alexiades Mathematics, University of Tennessee, Knoxville, TN 37996, USA and Oak Ridge National Laboratory, Oak Ridge TN 37831, USA E-mail address: [email protected] David Autrique Department of Chemistry, University of Antwerp, Antwerp, Belgium E-mail address: [email protected]