\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small Ninth MSU-UAB Conference on Differential Equations and Computational Simulations. \emph{Electronic Journal of Differential Equations}, Conference 20 (2013), pp. 93--101.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2013 Texas State University - San Marcos.} \vspace{9mm}} \begin{document}\setcounter{page}{93} \title[\hfilneg EJDE-2013/Conf/20/ \hfil Time-stepping for laser ablation] {Time-stepping for laser ablation} \author[H. Khanal, D. Autrique, V. Alexiades \hfil EJDE-2013/conf/20 \hfilneg] {Harihar Khanal, David Autrique, Vasilios Alexiades} % in alphabetical order \address{Harihar Khanal \newline Mathematics, Embry-Riddle Aeronautical University, Daytona Beach, FL 32114, USA} \email{Harihar.Khanal@erau.edu} \address{David Autrique \newline Physics, Optimas Research Center - TU Kaiserslautern, 67653 Kaiserslautern, Germany} \email{dautriqu@physik.uni-kl.de} \address{Vasilios Alexiades \newline Mathematics, University of Tennessee, Knoxville, TN 37996, USA} \email{alexiades@utk.edu} \thanks{Published October 31, 2013.} \subjclass[2000]{92C45, 35K60, 65M99} \keywords{Laser ablation; ode integration; time-stepper; peer methods; \hfill\break\indent plasma formation, collisional radiative model} \begin{abstract} Nanosecond laser ablation is a popular technique, applied in many areas of science and technology such as medicine, archaeology, chemistry, environmental and materials sciences. We outline a computational model for radiative and collisional processes occurring during ns-laser ablation, and compare the performance of various low and high order time-stepping algorithms. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{remark}[theorem]{Remark} \newtheorem{definition}[theorem]{Definition} \allowdisplaybreaks \section{Introduction}\label{S:1} The interaction of nanosecond pulsed lasers with solid targets and the properties of laser-produced plasmas have been investigated for many years. The laser energy deposited in a target can be used to remove sample material. This technique, known as laser ablation, is nowadays used in several medical, scientific and industrial applications, such as surgery, chemical analysis, pulsed laser deposition, laser remelting, laser machining, etc. \cite{Bauerle2011, Chrisey1994, Miller1998, Montaser1998}. In spite of the large number of scientific and practical applications of laser ablation, the mechanisms underlying laser-material interaction are not fully understood. In recent years, there have been significant strides in modelling and simulation of laser ablation using kinetic models \cite{kools1993monte, itina1997monte, garrelie1999study}, hydrodynamic models \cite{singh1990pulsed, aden1993laser, ho1995computational, le2000modeling, balazs1991expansion, anisimov2002selected, bulgakov1995dynamics, qaisar2003laser,wen2007expansion, mazhukin2007modeling, aghaei2008simulation, clair2011modelling, alexiades2010enthalpy, autrique2012multiphase}, as well as hybrid models such as \cite{harilal2003internal, itina2002laser}. We use a multiphase hydrodynamic model, described in the accompanying article \cite{autrique2013hydrodynamic}, which accounts for target heating, surface and volumetric mass removal, as well as plume expansion and plasma formation. The model consists of a tightly coupled system of partial differential equations (PDEs) and ordinary differential equations (ODEs). Due to the high computational cost, the problem needs efficient numerical algorithms. In this paper, which is a continuation of \cite{autrique2013hydrodynamic}, we present briefly a model accounting for the main collisional and radiative processes occurring during laser ablation, and compare the performance of various low- and high-order time-stepping algorithms. We implemented explicit Euler, non-adaptive RK of orders 2, 3, 4, adaptive RKF (RKFB4) and Dormand-Prince (DoPri5), as well as explicit PEER methods \cite{weiner2008explicit} (up to 9th order). After a brief description of the laser ablation process in \S \ref{S:2}, the mathematical model is outlined in \S\ref{S:3}. The computational approach is described and simulation results of various ODE time-stepping schemes (adaptive and non-adaptive methods of low and high orders) are presented in \S\ref{S:4}, and conclusions in \S\ref{S:5}. \section{Laser Ablation}\label{S:2} Laser ablation is characterized by several strongly coupled physical and chemical processes occurring in and above the target. The main processes can be summarized as follows. \subsection*{Target Heating} During the initial stage of ns-laser ablation, a part of the laser energy is absorbed in the vicinity of sample surface. The absorbed energy diffuses into the interior of the target, causing melting and removal (ablation) of some of the target material. Target heating is described by the simultaneous solution of an internal energy equation, continuity equation and a pressure relaxation equation, in a coordinate system attached to the ablating surface. Above the target, the evaporated particles achieve translational equilibrium within a few mean free paths by means of collisions, in a thin zone, known as the {\it Knudsen layer} (KL) \cite{Bauerle2011}. The Knudsen layer provides the connection between the target and the plasma plume. \subsection*{Plume Expansion and Vapor Flow} Beyond the Knudsen layer, the dense vapor plume ionizes during the laser action. A plasma is formed that absorbs laser energy, thus shielding the target surface from the incoming laser light. The absorption of laser energy by the plasma results in very high plume temperatures, velocities, species densities, pressures, etc. Afterwards, this hot plasma quickly expands into the ambient environment. The plasma plume is modeled by a set of compressible Euler equations \cite{autrique2013hydrodynamic}, closed by a multiphase equation of state (EOS) \cite{levashov2007eos}. This multiphase hydrodynamic model is presented in the accompanying article by Autrique et. al. \cite{autrique2013hydrodynamic}. As mentioned in \cite{autrique2013hydrodynamic}, species diffusion in the plume is not considered, due to the disparate time scales of the phenomena; diffusion processes become important much later, at near-s times \cite{le2000modeling} (here we simulate the process for 50 ns). \subsection*{Laser Induced Breakdown} In the irradiated vapor, several collisional and radiative processes take place. Optical breakdown is modeled by a dimensionless collisional radiative model. It involves a set of rate equations \cite{amoruso1999modeling, mazhukin2003optical, morel2010modeling} describing the temporal evolution of the electron density, electron and ion temperatures, as well as the atomic level populations in the plume. The rate equations account for single- and multi-photon ionization, heating due to inverse Bremsstrahlung absorption, radiative decay, electron impact excitation and ionization, as well as their respective recombination reactions \cite{chung2005flychk}. Following \cite{amoruso1999modeling, mazhukin2003optical, morel2010modeling}, the concentration of charged particles (ions and electrons) and the population of their excited levels is mathematically modeled by a (highly nonlinear) system of ordinary differential equations. The ODEs and associated reactions are listed in the Appendix \ref{S:appendix}. Combining equations \eqref{ode5}--\eqref{ode10}, the collisional radiative model constitutes an initial value problem for a system of ${\rm n}$ first order ODEs with m parameters, of the generic form: \begin{equation}\label{odesys} \frac{dy}{dt}= f(t,y,p(y)), \quad y(t_{\rm 0})=y_{\rm 0}\in \mathbb{R}^{\rm n}, \quad p\in \mathbb{R}^{\rm m}, \end{equation} where $y$ is an array containing the number densities of ions, electrons, excited states and their energies; $p$ is an array of parameters of the model, that contains all the rate constants (which depend on temperature, density and spectroscopic properties of the species). The ODE system \eqref{odesys} has to be solved simultaneously with the system of PDEs that accounts for target heating and vapor flow. An outline of the numerical procedures employed is given in \S \ref{S:4:1} below. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{ODE Integrators}\label{S:3} % We implemented various ODE integration schemes with fixed and with adaptive time-step sizes for the system of ODEs \eqref{odesys}. Due to the nature of the problem, and keeping future parallelization in mind, we confined our work to explicit methods only. In addition to explicit Euler (EE), we tested the low order Runge-Kutta methods RK2, RK3 and RK4. Then we tried the adaptive methods Runge-Kutta-Fehlberg (RKFB4) and Dormand-Prince (DoPri5). Finally, we implemented two-step {\it PEER methods} of orders 3 to 9, obtained from the EPPEER package \cite{Schmitt2012}. EPPEER \cite{Schmitt2012} is a Fortran 95 code for solving ODE initial value problems by explicit two-step peer methods with automatic step size control. %%%%%%%%%%%%%%%%%%%%%%%%%%% \subsubsection*{Explicit PEER Methods.} % At each time step of duration $h_{\rm k}$, from instant $t_{\rm k}$ to instant $t_{\rm {k+1}}=t_{\rm k}+h_{\rm k}$, a ${\rm s}$-stage PEER method solves an initial value problem % %% \[ \cred{y ^{\textbf{\text{\textquotesingle}}}} = f(t,y), \[ y ^\prime = f(t,y), \quad y(t_{\rm 0})=y_{\rm 0}, \] by constructing ${\rm s}$ approximations (stages) $Y_{\rm {k,i}}\approx y(t_{\rm {k,i}})$, ${\rm i=1,2},\dots {\rm s}$ at intermediate times $t_{\rm {k,i}}:=t_{\rm k}+h_{\rm k} c_{\rm i}$, with $t_{\rm {k,s}} = t_{\rm k}+h_{\rm k}$, given by \[ Y_{\rm {k,i}}=\sum_{\rm {j=1}}^{\rm s} b_{\rm {ij}}Y_{\rm {k-1,j}}+h_{\rm k} \sum_{\rm {j=1}}^{\rm s} \!\!a_{\rm {ij}}f(t_{\rm {k-1,j}},Y_{\rm {k-1,j}}) . \] The coefficients $a_{\rm {ij}}$, $b_{\rm {ij}}$, $c_{\rm i}$ are specific to each PEER method \cite{Schmitt2012}. At all stages ${\rm i=1,2},\dots {\rm s}$, the solutions $Y_{\rm {k,i}}$ possess the same accuracy and stability properties. \section{Numerical Simulations}\label{S:4} The system of ODEs \eqref{odesys} is solved using the time-stepping numerical schemes described in \S\ref{S:3}. This is done within the Finite Volume code for the system of PDEs describing target heating and plume expansion, see \cite{autrique2013hydrodynamic}. The overall algorithm is outlined below. \subsection{Computational Approach}\label{S:4:1} First, the heat conduction equation is solved in the target. The target and the plume domain are connected by jump relations that express the temperature, density and pressure variation across the Knudsen layer. The jump relations provide the inflow or outflow conditions for the Euler equations in the plume domain. As soon as the surface temperature exceeds the normal boiling point, material ends up in the plume domain and the ODE system \eqref{odesys} treating optical breakdown (plasma formation), is solved. Since one must resolve both radiative and collisional processes, the ODEs require a smaller time step than the PDEs. We implemented this as $\Delta t_{{\rm ODE}} = \Delta t_{\rm PDE} / factor$, with $factor$ an input parameter which we varied (see \S\ref{S:4:3}). As mentioned in \cite{autrique2013hydrodynamic}, after the breakdown stage, the collisional radiative model indicates that the (spatially averaged) plasma attains about equal electron, excitation, and heavy species temperatures. This means that a state close to Local Thermodynamic Equilibrium (LTE) is approached. At that instant, the rate equations are switched off and the temperature, electron density, and ion abundances are found from the Saha-Eggert equations \cite{Zeldovich2002}. \subsection{Simulation Setup}\label{S:4:2} We simulated laser ablation of copper (Cu) in a helium (He) background gas. The laser has a wavelength of 532~nm and a pulse width of 6~ns. Both target and background gas initially are in a stationary state at standard temperature and pressure. The ablation process is simulated for 50~ns. The entire setup (physical problem, parameter values, spatial grid, etc) and the code itself are identical with those used for the simulations reported in the accompanying article \cite{autrique2013hydrodynamic}, except here the peak intensity was set at {$I{\rm o} = 10^{13} \, {\rm W/m^2}$} for all runs. Only the time-stepper was changed from run to run, and the CPU timings were recorded. The numerical code is written in Fortran 90, compiled with Intel Fortran, and ran on Xeon-class processors (AMD Opteron 2378, 2400 MHz, 512 KB cache). \subsection{Simulation Results}\label{S:4:3} The comparison of CPU timings of various time-step\-pers employed in the code is presented in Table \ref{table1} and displayed in Figure \ref{fig1}. All time-steppers produce essentially identical results. The time step $factor$, mentioned in \S\ref{S:4:1}, is listed in the Table. Only the adaptive integrators RKFB4 and DoPri5 can use $factor$=10, all others require $factor$=20, i.e. smaller time-steps. We could not improve this factor by using higher order and adaptive time steppers, contrary to our expectations. Explicit Euler turns out to be somewhat faster than the other solvers. Moreover, fixed step size schemes perform better than the adaptive ones. \begin{table}[ht] \caption{CPU timings of ablation processes for 50~ns} \begin{tabular}{|l|c|c|c|c|c|} \hline Time-stepper & Order & $\Delta$t $factor$ & Step size & CPU (sec)\\ \hline \hline EE & 1 & 20 & fixed & 7288 \\ \hline RK2 & 2 & 20 & fixed & 7353 \\ \hline RK3 & 3 & 20 & fixed & 7398 \\ \hline RK4 & 4 & 20 & fixed & 7453 \\ \hline RKFB4 & 4 & 10 & adaptive & 7412 \\ \hline DoPri5 & 5 & 10 & adaptive & 7442 \\ \hline Peer & 3 & 20 & adaptive & 7949 \\ \hline Peer & 4 & 20 & adaptive & 8442 \\ \hline Peer & 5 & 20 & adaptive & 8468 \\ \hline Peer & 6 & 20 & adaptive & 9002 \\ \hline Peer & 7 & 20 & adaptive & 9026 \\ \hline Peer & 8 & 20 & adaptive & 9519 \\ \hline Peer & 9 & 20 & adaptive & 9559 \\ \hline \end{tabular} \label{table1} \end{table} \begin{figure}[ht] \begin{center} \includegraphics[width=0.6\textwidth]{fig1} % fig_timings.pdf \end{center} \caption{CPU timings of laser ablation for 50~ns using various ODE integrators (ee: Explicit Euler, rkN: Runge-Kutta of order N, rkf: Runge-Kutta-Fehlberg, dp: Dormand-Prince, pN: peer method of order N)} \label{fig1} \end{figure} \section{Conclusions}\label{S:5} We outlined a collisional radiative model for ns-laser induced breakdown in an expanding copper plume. The model accounts for photon ionization, electron impact excitation and ionization, as well as the respective recombination reactions. In a next step, numerical simulations were performed for ns-laser ablation of a copper target at a laser wavelength of 532~nm and peak intensity $10^{13} \, {\rm W/m^2}$, comparing the performance of 13 ODE integrators. We found that lower order methods are faster than higher order methods and Explicit Euler (EE) turns out to be the fastest. This may be attributed to the very high cost of evaluating the right-hand-sides of the ODEs; thus, the fewer the evaluations, the better the performance. Surprisingly, the higher order accuracy of high order schemes did not translate to speedup, as we had to use the fine step size to maintain accuracy. Finally, contrary to our expectations, non-adaptive methods turned out to perform better than adaptive ones, and the PEER methods are slower than the other integrators. Further tests are under way, exploring other solvers and coding strategies. \section{Appendix}\label{S:appendix} The following relations hold for both vapor (copper) and background gas (helium). The main radiative and collisional processes, as well as their corresponding ODEs are listed below, compiled from \cite{chung2005flychk}, \cite{mazhukin2003optical}, \cite{morel2010modeling}. The number of energy levels in oxidation (charge) state ${\rm k}$ is denoted by ${\rm n(k)}$, and the running indices ${\rm i}$ and ${\rm j}$ denote energy levels. $N_{\rm {j}}^{\rm {k}}$ is the number density of species of charge state ${\rm k}$ excited to level ${\rm j}$, $N^{{\rm k}}$ is the total ion number density in charge state ${\rm k}$, and $N_e$ the electron number density. The electrons that participate in the various processes can be distinguished by their energy $\epsilon $. $\Delta {{E}_{{\rm ij}}}$ is the energy difference between levels ${\rm i}$ and ${\rm j}$, where the indices ${\rm i}$ and ${\rm j}$ can belong to the same oxidation state ${\rm k}$, or where ${\rm i}$ can belong to state ${\rm k}$ and ${\rm j}$ to ${\rm k+1}$, respectively. The corresponding angular frequencies are denoted as $\nu_{{\rm ij}}$ ( = $\Delta {{E}_{\rm {ij}}} / \hbar$ ), whereas $\omega_{{\rm las}}$ designates the angular laser frequency. The various rate coefficients and electron energy source terms are denoted by $k_{\rm {LAB,i,j}}$ and $S_{\rm {LAB,j}}^{\rm {k}}$ for various labels \lq\lq ${\rm LAB}$\rq\rq, respectively: ${\rm PI}$ = photo ionization, ${\rm dec}$ = radiative decay, ${\rm exc}$ = excitation, ${\rm dexc}$ = de-excitation, ${\rm ei}$ = electron impact ionization, ${\rm tbod}$ = three-body recombination. Other symbols appearing below are: $\hbar$ = Dirac constant, $m_{\rm e}$ = mass of electron, $m_{\rm h}$ = mass of heavy species, $m$ = number of photons needed for a certain transition ${\rm i}$ of state ${\rm k}$ to ${\rm j}$ of state ${\rm k+1}$. \subsection*{I. Radiative Processes} \quad \\ $\bullet$ Laser-induced photo ionization \& recombination: ${\rm Cu_i^k}+ m \hbar\omega{\rm las} \leftrightarrow {\rm {Cu_0}^{k+1}+ e(\epsilon) }$ \begin{equation}\label{ode1} \begin{gathered} \frac{dN_{\rm {j}}^{\rm {k}}}{dt} = \sum_{\rm {i=1}}^{\rm {n(k-1)}}{k_{\rm {PI,i,j}}^{\rm {k-1}}N_{\rm {i}} ^{\rm {k-1}}}-\sum_{\rm {i=1}}^{\rm {n(k+1)}}{k_{\rm {PI,i,j}}^{\rm {k+1}} N_{\rm {j}}^{\rm {k}}} \,, \\ S_{\rm {PI,j}}^{\rm {k}} = \sum_{\rm {i=1}}^{\rm {n(k-1)}}{k_{\rm {PI,i,j}}^{\rm {k-1}}N_{\rm {i}} ^{\rm {k-1}}(m \hbar{\omega{\rm las}}-\Delta {{E}_{\rm {ij}}})} +\sum_{\rm {i=1}}^{\rm {n(k+1)}}{k_{\rm {PI,i,j}}^{\rm {k+1}} N_{\rm {j}}^{\rm {k}}}(m \hbar{\omega{\rm las}}-\Delta {{E}_{\rm {ij}}}). \end{gathered} \end{equation} $\bullet$ Radiative decay: $ {\rm Cu_j^k} \to {\rm Cu_i^{k}}+\hbar\omega_{\rm ij} $ \begin{equation} \frac{dN_{\rm {j}}^{\rm {k}}}{dt}=\sum_{\rm {i=j+1}}^{\rm {n(k)}}{k_{\rm {dec,i,j}}^{\rm {k}}N_{\rm {i}}^{\rm {k}}}-\sum_{\rm {i=1}}^{\rm {j-1}}{k_{\rm {dec,i,j}}^{\rm {k}}N_{\rm {j}}^{\rm {k}}} \, . \label{ode2} \end{equation} \subsection*{II. Collisional Processes} \quad\\ $\bullet$ Collisional excitation \& de-excitation: $ {\rm Cu_i^k+ e(\epsilon)} \leftrightarrow {\rm Cu_j^{k} + e(\epsilon') }$ \begin{equation} \label{ode3} \begin{aligned} \frac{dN_{\rm {j}}^{\rm {k}}}{dt} &=\sum_{\rm {i=1}}^{\rm {j-1}} {\Big( k_{\rm {exc,i,j}}^{\rm {k}}N_{\rm {i}}^{\rm {k}}-k_{\rm {dexc,i,j}} ^{\rm {k}}N_{\rm {j}}^{\rm {k}} \Big)}{{N}_{\rm {e}}}\\ &\quad +\sum_{\rm {i=j+1}}^{\rm {n(k)}}{\Big( k_{\rm {dexc,i,j}}^{\rm {k}}N_{\rm {i}}^{\rm {k}}-k_{\rm {exc,i,j}}^{\rm {k}} N_{\rm {j}}^{\rm {k}} \Big)}{{N}_{\rm {e}}} \, , \\ S_{\rm {ex,j}}^{\rm {k}} &=-\sum_{\rm {i=1}}^{\rm {j-1}} {\Big( k_{\rm {exc,i,j}}^{\rm {k}}N_{\rm {i}}^{\rm {k}} -k_{\rm {dexc,i,j}}^{\rm {k}}N_{\rm {j}}^{\rm {k}} \Big)} \Delta {{E}_{\rm {ij}}}{{N}_{\rm {e}}}\\ &\quad +\sum_{\rm {i=j+1}}^{\rm {n(k)}} {\Big(k_{\rm {dexc,i,j}}^{\rm {k}}N_{\rm {i}}^{\rm {k}} -k_{\rm {exc,i,j}}^{\rm {k}}N_{\rm {j}}^{\rm {k}} \Big)} \Delta {{E}_{\rm {ij}}}{{N}_{\rm {e}}} \, . \end{aligned} \end{equation} $\bullet$ Electron impact ionization \& three-body recombination: \\ $ {\rm Cu_i^{k}+ e(\epsilon)} \leftrightarrow {\rm Cu_j^{k+1} + e(\epsilon') +e(\epsilon'')}$ \begin{equation} \label{ode4} \begin{aligned} \frac{dN_{\rm {j}}^{\rm {k}}}{dt} &=\sum_{\rm {i=1}}^{\rm {n(k-1)}} {\Big( k_{\rm {ei,i,j}}^{\rm {k-1}}N_{\rm {i}}^{\rm {k-1}} -k_{\rm {tbod,i,j}}^{\rm {k-1}}N_{\rm {j}}^{\rm {k}}{{N}_{\rm {e}}} \Big)} {{N}_{\rm {e}}}\\ &\quad +\sum_{\rm {i=1}}^{\rm {n(k+1)}}{\Big( k_{\rm {tbod,i,j}} ^{\rm {k+1}}{{N}_{\rm {e}}}N_{\rm {i}}^{\rm {k+1}} -k_{\rm {ei,i,j}}^{\rm {k+1}}N_{\rm {j}}^{\rm {k}} \Big)}{{N}_{\rm {e}}} \, , \\ S_{\rm {ei,j}}^{\rm {k}} &=-\sum_{\rm {i=1}}^{\rm {n(k-1)}} {\Big( k_{\rm {ei,i,j}}^{\rm {k-1}}N_{\rm {i}}^{\rm {k-1}} -k_{\rm {tbod,i,j}}^{\rm {k-1}}N_{\rm {j}}^{\rm {k}}{{N}_{\rm {e}}} \Big)} \Delta {{E}_{\rm {ij}}}{{N}_{\rm {e}}}\\ &\quad +\sum_{\rm {i=1}}^{\rm {n(k+1)}} {\Big( k_{\rm {tbod,i,j}}^{\rm {k+1}}{{N}_{\rm {e}}}N_{\rm {i}}^{\rm {k+1}} -k_{\rm {ei,i,j}}^{\rm {k+1}}N_{\rm {j}}^{\rm {k}} \Big)} \Delta {{E}_{\rm {ij}}}{{N}_{\rm {e}}} \, . \end{aligned} \end{equation} \subsection*{III. The ODE System} The ODEs characterizing the overall evolution of concentration of species, and of temperatures of electrons and ions are presented below. Here $m_{\rm e}$ and $m_{\rm h}$ , $T_{\rm e}$ and $T_{\rm h}$, ${N}_{\rm e}$ and ${N}_{\rm h}$ denote the atomic mass, temperature and total density of the electrons and heavy species, respectively. The electron-ion and electron-neutral collision frequencies are denoted by $\nu_{\rm {ei}}$ and $\nu_{\rm {en}}$, respectively. ${S}_{\rm {IB}}$ is the electron energy source term due to inverse Bremsstrahlung absorption. The system is solved up to a maximum charge state of 2. \begin{equation} \label{ode5} \begin{aligned} \frac{dN_{\rm {j}}^{\rm {k}}}{dt} &=\sum_{\rm {i=1}}^{\rm {j-1}} {\left( k_{\rm {exc,i,j}}^{\rm {k}}N_{\rm {i}}^{\rm {k}}-k_{\rm {dexc,i,j}} ^{\rm {k}}N_{\rm {j}}^{\rm {k}} \right)}{{N}_{\rm {e}}} + \sum_{\rm {i=j+1}}^{\rm {n(k)}} {\left( k_{\rm {dexc,i,j}}^{\rm {k}}N_{\rm {i}}^{\rm {k}} -k_{\rm {exc,i,j}}^{\rm {k}}N_{\rm {j}}^{\rm {k}} \right)} {{N}_{\rm {e}}} \\ &\quad - \sum_{\rm {i=1}}^{\rm {j-1}}{k_{\rm {dec,i,j}}^{\rm {k}}N_{\rm {j}} ^{\rm {k}}}+\sum_{\rm {i=j+1}}^{\rm {n(k)}}{k_{\rm {dec,i,j}}^{\rm {k}} N_{\rm {i}}^{\rm {k}} } \\ &\quad + {\sum_{\rm {i=1}}^{\rm {n(k-1)}}{\left( k_{\rm {ei,i,j}} ^{\rm {k-1}}N_{\rm {i}}^{\rm {k-1}}-k_{\rm {tbod,i,j}} ^{\rm {k-1}}N_{\rm {j}}^{\rm {k}}{{N}_{\rm {e}}} \right)} {{N}_{\rm {e}}}} \\ &\quad + \sum_{\rm {i=1}}^{\rm {n(k+1)}} {\left( k_{\rm {tbod,i,j}}^{\rm {k+1}}{{N}_{\rm {e}}} N_{\rm {i}}^{\rm {k+1}}-k_{\rm {ei,i,j}}^{\rm {k+1}}N_{\rm {j}} ^{\rm {k}} \right)}{{N}_{\rm {e}}} \\ &\quad +\sum_{\rm {i=1}}^{\rm {n(k-1)}}{k_{\rm {PI,i,j}} ^{\rm {k-1}}N_{\rm {i}} ^{\rm {k-1}}-} \sum_{\rm {i=1}}^{\rm {n(k+1)}} {k_{\rm {PI,i,j}} ^{\rm {k+1}}N_{\rm {j}}^{\rm {k}}} \,, \end{aligned} \end{equation} \begin{gather} {N}^{\rm {k}}=\sum_{\rm {j=1}}^{\rm {n(k)}}{N_{\rm {j}}^{\rm {k}}} \quad, {N}_{\rm {h}}=\sum_{\rm {k=0}}^{\rm {2}}{{{N}^{\rm {k}}}} \quad, {N}_{\rm {e}}=\sum_{\rm {k=0}}^{\rm {2}}{k{{N}^{\rm {k}}}} \, , \label{ode8} \\ \frac{d(\frac{3}{2}kT_{\rm {h}}N_{\rm h})}{dt} =3\frac{{{m}_{\rm {e}}}}{{{m}_{\rm {h}}}} \left( {{T}_{\rm {e}}}-{{T}_{\rm {h}}} \right)\left( {{\nu}_{\rm {en}}} +{{\nu}_{\rm {ei}}} \right){{N}_{\rm {e}}} \, , \label{ode9} \\ \frac{d(\frac{3}{2}kT_{\rm e}N_{\rm e})}{dt} =-\frac{d}{dt}\left(\! \tfrac{3}{2}k{{T}_{\rm {h}}}{{N}_{\rm {h}}} \right) +{{S}_{\rm {IB}}} + \sum_{\rm {k=0}}^{\rm {2}} {\sum_{\rm {j=1}}^{\rm {n(k)}}{\left( S_{\rm {ei,j}}^{\rm {k}} + S_{\rm {ex,j}}^{\rm {k}}+S_{\rm {PI,j}}^{\rm {k}} \right)}} \, . \label{ode10} \end{gather} \subsection*{Acknowledgments} The authors wish to thank O. Rosmej for her valuable comments and efforts. P. Levashov, K. Khishenko and M. Povarnitsyn are acknowledged for the equation-of-state data set and advice. The second author acknowledges financial support from the Deutsche Forschungsgemeinschaft (Emmy Noether-\-Program, grant RE 1141/11). \begin{thebibliography}{10} \bibitem{aden1993laser} M.~Aden, E. W. Kreutz, A.~Voss; \emph{Laser-induced plasma formation during pulsed laser deposition}, Journal of Physics D: Applied Physics \textbf{26} (1993), 1545--1553. \bibitem{aghaei2008simulation} M.~Aghaei, S.~Mehrabian, S. H. Tavassoli; \emph{Simulation of nanosecond pulsed laser ablation of copper samples: A focus on laser induced plasma radiation}, Journal of Applied Physics \textbf{104} (2008), no.~5, 053303. \bibitem{alexiades2010enthalpy} V.~Alexiades, D.~Autrique; \emph{Enthalpy model for heating, melting, and vaporization in laser ablation}, Electronic Journal of Differential Equations \textbf{Conf. 19} (2010), 1--14. \bibitem{amoruso1999modeling} S.~Amoruso; \emph{Modeling of {UV} pulsed-laser ablation of metallic targets}, Applied Physics A: Materials Science \& Processing \textbf{69} (1999), no.~3, 323--332. \bibitem{anisimov2002selected} S. I. Anisimov, B. S. Luk'yanchuk; \emph{Selected problems of laser ablation theory}, Physics-Uspekhi \textbf{45} (2002), no.~3, 293--324. \bibitem{autrique2013hydrodynamic} D.~Autrique, V.~Alexiades, H.~Khanal; \emph{Hydrodynamic modeling of ns-laser ablation}, Electronic Journal of Differential equations, conference 20 (2013), pp. 1--14. \bibitem{autrique2012multiphase} D.~Autrique, Z.~Chen, V.~Alexiades, A.~Bogaerts, B.~Rethfeld; \emph{A multiphase model for pulsed ns-laser ablation of copper in an ambient gas}, AIP Conference Proceedings, vol. 1464, 2012, p.~648. \bibitem{balazs1991expansion} L.~Balazs, R.~Gijbels, and A.~Vertes; \emph{Expansion of laser-generated plumes near the plasma ignition threshold}, Analytical chemistry \textbf{63} (1991), no.~4, 314--320. \bibitem{Bauerle2011} D.~B\"{a}uerle; \emph{{L}aser {P}rocessing and {C}hemistry}, Springer Verlag, Berlin, 2011. \bibitem{bulgakov1995dynamics} A. V. Bulgakov and N. M. Bulgakova; \emph{Dynamics of laser-induced plume expansion into an ambient gas during film deposition}, Journal of physics D: applied physics \textbf{28} (1995), no.~8, 1710--1718. \bibitem{Chrisey1994} D. B. Chrisey, G. K. Hubler; \emph{{P}ulsed {L}aser {D}eposition of {T}hin {F}ilms}, Wiley, New York, 1994. \bibitem{chung2005flychk} H. K. Chung, M. H. Chen, W. L. Morgan, Y.~Ralchenko, R. W. Lee; \emph{{FLYCHK}: {G}eneralized population kinetics and spectral model for rapid spectroscopic analysis for all elements}, High Energy Density Physics \textbf{1} (2005), no.~1, 3--12. \bibitem{clair2011modelling} G.~Clair, D.~L'Hermite; \emph{{1D modelling of nanosecond laser ablation of copper samples in argon at P = 1 atm with a wavelength of 532 nm}}, Journal of Applied Physics \textbf{110} (2011), no.~8, 083307. \bibitem{garrelie1999study} F.~Garrelie, C.~Champeaux, A.~Catherinot; \emph{Study by a {M}onte {C}arlo simulation of the influence of a background gas on the expansion dynamics of a laser-induced plasma plume}, Applied Physics A: Materials Science \& Processing \textbf{69} (1999), no.~1, 45--50. \bibitem{harilal2003internal} S. S. Harilal, C. V. Bindhu, M. S. Tillack, F.~Najmabadi, A.C. Gaeris; \emph{Internal structure and expansion dynamics of laser ablation plumes into ambient gases}, Journal of Applied Physics \textbf{93} (2003), no.~5, 2380--2388. \bibitem{ho1995computational} J. R. Ho, C. P. Grigoropoulos, J. A. C. Humphrey; \emph{Computational study of heat transfer and gas dynamics in the pulsed laser evaporation of metals}, Journal of Applied Physics \textbf{78} (1995), no.~7, 4696--4709. \bibitem{itina2002laser} T. E. Itina, J.~Hermann, P.~Delaporte, M.~Sentis; \emph{Laser-generated plasma plume expansion: {C}ombined continuous-microscopic modeling}, Physical Review E \textbf{66} (2002), no.~6, 066406. \bibitem{itina1997monte} T. E. Itina, V. N. Tokarev, W.~Marine, M.~Autric; \emph{Monte {C}arlo simulation study of the effects of nonequilibrium chemical reactions during pulsed laser desorption}, Journal of Chemical Physics \textbf{106} (1997), no.~21, 8905--8912. \bibitem{kools1993monte} J. C. S. Kools; \emph{Monte {C}arlo simulations of the transport of laser-ablated atoms in a diluted gas}, Journal of Applied Physics \textbf{74} (1993), no.~10, 6401--6406. \bibitem{le2000modeling} H. C. Le, D. E. Zeitoun, J. D. Parisse, M.~Sentis, W.~Marine; \emph{Modeling of gas dynamics for a laser-generated plasma: {P}ropagation into low-pressure gases}, Physical Review E \textbf{62} (2000), no.~3 Pt B, 4152--4161. \bibitem{levashov2007eos} P. R. Levashov, K. V. Khishchenko; \emph{{ITTEOS} 5.8 software for calculation of {EOS} for {M}etals}, 2007. \bibitem{mazhukin2003optical} V. I. Mazhukin, V. V. Nossov, M. G. Nickiforov, I.~Smurov; \emph{Optical breakdown in aluminum vapor induced by ultraviolet laser radiation}, Journal of Applied Physics \textbf{93} (2003), no.~1, 56--66. \bibitem{mazhukin2007modeling} V. I. Mazhukin, V. V. Nossov, I.~Smurov; \emph{Modeling of plasma-controlled evaporation and surface condensation of {A}l induced by 1.06 and 0.248 µm laser radiations}, Journal of Applied Physics \textbf{101} (2007), no.~2, 024922--024922. \bibitem{Miller1998} J. C. Miller, R. F. Haglund; \emph{Laser ablation and desorption}, Academic Press, New York, 1998. \bibitem{Montaser1998} A.~Montaser; \emph{{I}nductively {C}oupled {P}lasma {M}ass {S}pectrometry}, Wiley, New York, 1998. \bibitem{morel2010modeling} V.~Morel, A.~Bultel, B. G. Ch{\'e}ron; \emph{Modeling of thermal and chemical non-equilibrium in a laser-induced aluminum plasma by means of a {C}ollisional-{R}adiative model}, Spectrochimica Acta, Part B: Atomic Spectroscopy \textbf{65} (2010), no.~9, 830--841. \bibitem{qaisar2003laser} M. S. Qaisar, G. J. Pert; \emph{Laser ablation of {M}g, {C}u, and {P}b using infrared and ultraviolet low-fluence lasers}, Journal of Applied Physics \textbf{94} (2003), no.~3, 1468--1477. \bibitem{Schmitt2012} B. A. Schmitt, R.~Weiner; \emph{Manual for the {E}xplicit {P}arallel {P}eer {C}ode {EPPEER}}, \\ {www.mathematik.uni-marburg.de/~schmitt/peer}, 2012. \bibitem{singh1990pulsed} R. K. Singh, J.~Narayan; \emph{Pulsed-laser evaporation technique for deposition of thin films: {P}hysics and theoretical model}, Physical Review B \textbf{41} (1990), no.~13, 8843--8859. \bibitem{weiner2008explicit} R.~Weiner, K.~Biermann, B. A. Schmitt, H.~Podhaisky; \emph{Explicit two-step peer methods}, Computers \& Mathematics with Applications \textbf{55} (2008), no.~4, 609--619. \bibitem{wen2007expansion} S. B. Wen, X.~Mao, R.~Greif, R. E. Russo; \emph{Expansion of the laser ablation vapor plume into a background gas. {I}. {A}nalysis}, Journal of Applied Physics \textbf{101} (2007), no.~2, 023114. \bibitem{Zeldovich2002} B.~Zel'dovich, Y. P. Raizer; \emph{Physics of {S}hock {W}aves and {H}igh-{T}emperature {H}ydrodynamic {P}henomena}, vol. 1,2, Dover, New York, 2002. \end{thebibliography} \end{document}