\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}