\documentclass[reqno]{amsart}
%\usepackage[notref,notcite]{showkeys}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
Eighth Mississippi State - UAB Conference on Differential Equations and
Computational Simulations.
{\em Electronic Journal of Differential Equations},
Conference 19 (2010),  pp. 1--14.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2010 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document} \setcounter{page}{1}
\title[\hfilneg EJDE-2010/Conf/19/\hfil Enthalpy model]
{Enthalpy model for heating, melting, and vaporization in laser ablation}

\author[V. Alexiades, D. Autrique \hfil EJDE/Conf/19 \hfilneg]
{Vasilios Alexiades, David Autrique} 

\address{Vasilios Alexiades \newline
Mathematics, University of Tennessee,
 Knoxville, TN 37996, USA \newline
 \emph{and}  Oak Ridge National Laboratory, Oak Ridge TN 37831, USA}
\email{alexiades@utk.edu}

\address{David Autrique \newline
Department of Chemistry, University of Antwerp,
 Antwerp, Belgium}
 \email{david.autrique@ua.ac.be}

\thanks{Published September 25, 2010.}
\subjclass[2000]{92C45, 35K60, 65M99}
\keywords{Laser ablation; phase change; enthalpy scheme; thermochemistry; 
\hfill\break\indent
finite volume scheme}

\begin{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. 
\end{abstract}

\maketitle
\numberwithin{equation}{section}

\section{Introduction}\label{S:1}


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 
\cite{ready71,vonAllmen87,chrisey94}.

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 
\cite{vertes93,broekaert02,becker02,gunther02}.

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. \ref{fig:1}.

\begin{figure}[ht]
 \centering
 \includegraphics[width=0.95\textwidth]{fig1}  %  height=2.3in
\label{fig:1} \caption{
   Schematic of the Laser Ablation - Inductively Coupled Plasma 
   Mass Spectrometry (LA-ICP-MS) process. (Courtesy of C.C. Garcia, 
   ISAS, Dortmund)}
\end{figure}

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
\cite{noorbatcha87,ellegaard99,itina97, garrelie98,kools93},
hydrodynamic
\cite{bogaerts03,chen05,gusarov00,bleiner05,lindner09},
or hybrid models \cite{itina02,itina03}.
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 \S\ref{S: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 \S\ref{S:2} we briefly describe some issues that arise during the 
modeling of \emph{nanosecond laser ablation}.  
In \S\ref{S: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 \S\ref{S:4},
and final conclusions in \S\ref{S:5}.


\section{Nanosecond laser ablation}\label{S:2}

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 \emph{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
\cite{vonAllmen87,bleiner06}.
An effective hydrodynamic model of the entire ablation process 
should posess the following highly desirable features:

\begin{itemize}
 \item seamless treatment of heating / melting / vaporization,
 \item temperature-dependent thermophysical and optical properties,
 \item treatment of vapor as a non-ideal gas by means of an equation of 
	state(EOS) up to and above the critical temperature $T_{\rm crit}~$,
 \item treatment of the plasma-sheet at the target surface,
 \item coupling target and plume by means of the Knudsen layer, 
 \item treatment of plume expansion and recondensation, in vacuum or in a background gas
 \item dealing with plasma formation (ions and electrons) mechanisms,
 \item treating absorption of laser energy in the plume and shielding of the target,
 \item capturing (very) strong shocks via high resolution numerical schemes,
 \item be applicable in 1D, 2D axisymmetric and 3D spatial dimensions.
\end{itemize}

\noindent
Additional challenges in modeling laser ablation arise from:

\begin{itemize}
 \item extreme space and time scales,
 \item extreme gradients: temperature may rise to thousands of degrees locally,
 \item extreme variation in thermophysical properties,
 \item 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 $T_{\rm crit}$ ( $\sim 8000$ K ),
 \item the need for T-dependent and wavelength dependent  optical data.
\end{itemize}

\section{Enthalpy model for target heating - melting - vaporization}\label{S:3}

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$.
\begin{equation}
 \partial{}_t H + \nabla \cdot \vec{Q} = S \, .
\end{equation}\label{eq:3:1}
%
The laser source term $S$, at time $t$ and position $z$, is given by
%
\begin{equation}
  S(t, z) = I(t) (1-R) \, \alpha e^{- \alpha z}  ~,
\end{equation}\label{eq:3:2}
%
with $R$ and $\alpha$ the reflection and absorption coefficients of the 
target, respectively, and intensity
%
\begin{equation}
  I(t) = \beta I_o \, e^{ - (t-t_{max})^2 / 2\sigma^2 }\, .
\end{equation}\label{eq:3:3}
%
Here $\beta$ 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 \S\ref{S: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} {\mbox W/m^2}$, 
at a wavelength of 266 nm. 
In Fig. \ref{fig:2} the laser intensities before ($\beta = 1$) and 
after shielding at the surface are shown.
% 
 \begin{figure}[ht]	
 \centering
 \includegraphics[width=0.95\textwidth]{fig2} %,height=4.0in
 \label{fig:2}\caption{
	Time course $I(t)$ of laser pulse intensity with FWHM=10 ns, 
	peak $7 \times 10^{12} \mbox W/m^2$.
     }
\end{figure}

There are three time scales that characterize laser heating: 
$\tau_e = C_e / \gamma$, the electron cooling time, 
$\tau_i=C_i / \gamma$, the lattice heating time, and
$\tau_l$, the duration of laser pulse. 
Here $C_e$ and $C_i$ are the heat capacities of the lattice and 
the electron subsystems, respectively, and
$\gamma$ a parameter characterizing the electron-lattice coupling in 
the target. Since $\tau_e$ is of the order of 1 ps and 
$\tau_l$ is of the order of ns, the electron absorbed energy has 
enough time to be transferred to the lattice 
\cite{chichkov96}.
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 
\cite{harrington67}.
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, $\vec{Q} = - k \nabla 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 $T_m$ as well as above the normal 
boiling point $T_b$. 
Therefore the vaporization temperature $T_V (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 \eqref{eq:3:1} 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.

\subsection{Phases from volumetric enthalpy}\label{S:3:1}

The development is based on the basic Gibbs relation for the exact
differential of enthalpy \cite{book} in any phase $i = S, L, V$ ~:
%
\begin{equation}
 d H^i = {C_p}^i dT + [ -{\alpha}^i T + 1 ] dP \quad i=S, L, V~,
\end{equation}\label{eq:3:4}
%
where 
$C_p^i$ and $\alpha^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 \eqref{eq:3:4}, relative to a \emph{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 "${}_{\mbox o}$" 
denote interfacial (switch) values.  
Clearly, $H = \rho 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
\emph{volume} phase fraction, $\Lambda^i$, 
or by its \emph{mass} phase fraction, $\lambda^i$. 
The fractions are related by: $\rho^i \Lambda^i = \rho \lambda^i$
($i = S,L,V$), where $\rho$ = total density (of mixture) = 
$\Sigma_i \rho^i \Lambda^i$.  Note that fractions must add up to 1: 
$\Lambda^S + \Lambda^L + \Lambda^V = 1$ and 
$\lambda^S + \lambda^L + \lambda^V = 1$.
%
For simplicity, we assume the density of solid, $\rho^S$, and of 
liquid, $\rho^L$, remain constant (but $\rho^S \ne \rho^L$). 

As reference state we choose \emph{solid} at temperature $T_m$ and
pressure $P_{\rm ref}=1$ atm. The value $h_{\rm ref}$ 
(and thus also $H_{\rm ref} = \rho^S h_{\rm ref}$) of the enthalpy
at the reference state can be arbitrary; it will be chosen later
(so as to make $\Delta H^{vap} > 0$).
By integration of \eqref{eq:3:4} within each phase, and changing phase 
at ($T_m , P_{\rm ref}$) (from solid to liquid)
and at ($T_V , P_{\rm ref}$) (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 {\bf\emph{switch}} values
$H_o^S$, $H_o^L$, $H_o^{LV}(T_V)$, and $H_o^V (T_V)$, 
defined below, with
$T_V$ any desired vaporization temperature ($T_m < T_V \leq T_b$).

\vspace{2mm}
\noindent
$\bullet$ \textbf{Solid phase}\quad ($T \leq T_m$): \quad
$H < {\bf H_o^S} := H_{\rm ref} \equiv \rho^S h_o^S$, (set $h_o^S := h_{\rm ref}$).
\begin{center} 
Then $\Lambda^S = 1$, $\Lambda^L = 0$, $\Lambda^V = 0$, 
$\rho = \rho^S$, \quad
$H(T) = H_o^S + \int_{T_m}^T {C_p^{S} (\tau)} d\tau $. \\ 
Thus, given $H$, we can find the temperature $T$ ($< T_m$) by solving 
the equation:\\
$\int_{T_m}^T {C_p^{S} (\tau)} d\tau = H - H_o^S$ \quad for $T$.
\end{center} 

\noindent
$\bullet$ \textbf{Solid-Liquid transition}\quad (at $T = T_m$):\quad
${\bf H_o^S} < H < {\bf H_o^L} := \rho^L h_o^L$, \\
with $h_o^L = h_o^S + \Delta h^{fus}$.  Set 
$\Delta H^{fus} := H_o^L - H_o^S \equiv (\rho^L-\rho^S ) h_o^S + \rho^L \Delta h^{fus}$,
which can be made $> 0$ by choosing $h_{\rm ref}$ to be $\leq 0$.
%
\begin{center} 
Then $T = T_m$, \quad $\Lambda^L = (H-H_o^L)/\Delta H^{fus}$, 
$\Lambda^S = 1-\Lambda^L$, $\Lambda^V = 0$, 
$\rho = \Lambda^S \rho^S + \Lambda^L \rho^L$,
$H = H_o^L + \Lambda^L \Delta H^{fus}$, $h = H/\rho$.
\end{center} 

\noindent
$\bullet$ \textbf{Liquid phase}\quad ($T_m < T < T_V$):
${\bf H_o^L} < H < {\bf H_o^{LV}(T_V)} := H_o^L 
+ \int_{T_m}^{T_V} {C_p^{L} (\tau)} d\tau $.
%
\begin{center} 
Then $\Lambda^L = 1$, $\Lambda^S = 0$, $\Lambda^V = 0$, 
$\rho = \rho^L$,
$H(T) = H_o^L + \int_{T_m}^T {C_p^{L} (\tau)} d\tau $. \\
Thus, given $H$, we can find the temperature $T$ ($T_m < T < T_V$) 
by solving for $T$ the equation:\quad 
$\int_{T_m}^T {C_p^{L} (\tau)} d\tau = H - H_o^L$.
\end{center} 

\noindent
$\bullet$ \textbf{Liquid-Vapor transition}\quad (at $T = T_V$):\quad
${\bf H_o^{LV}(T_V)} < H < {\bf H_o^V (T_V, P_{\rm ref}) } :=\rho^V(T_V, P_{\rm ref}) h_o^V$,
with $h_o^V (T_V, P_{\rm ref}) := h_o^{LV}(T_V) 
+ \Delta h^{vap}(T_V, P_{\rm ref})$, where \\
$\rho^V(T_V, P_{\rm ref})$ = vapor density at ($T_V, P_{\rm ref}$), 
and $\Delta h^{vap}(T_V, P_{\rm ref})$ = heat of vaporization at 
($T_V, P_{\rm ref}$). Set 
$\Delta H^{vap} := H_o^V - H_o^{LV} \equiv (\rho^V-\rho^L ) h_o^{LV}(T_V) + \rho^V \Delta h^{vap}(T_V, P_{\rm ref})$, 
which can be made $> 0$  for any $T_m < T < T_{\rm crit}$ 
by choosing $h_{\rm ref}$ sufficiently negative, namely 
$h_{\rm ref} < -\Delta h^{fus} - c_p^L \, [T_{\rm crit} - T_m]$, 
thanks to the fact that $\Delta h^{vap}(T_{\rm crit}) = 0$.
%
\begin{center} 
Then $\Lambda^V = (H-H_o^{LV})/\Delta H^{vap}$, $\Lambda^L = 1-\Lambda^V$, $\Lambda^S = 0$, 
$\rho = \Lambda^L \rho^L + \Lambda^V \rho^V$, 
$H = \rho h = H_o^{LV}(T_V) + \Lambda^V \Delta H^{vap}$,
$h=H/\rho$. \\ 
\end{center} 

\noindent
$\bullet$ \textbf{Vapor phase}\quad ($T > T_V$):\quad
$H > {\bf H_o^V (T_V, P_{\rm ref}) }$.
\begin{center} 
Then $\Lambda^V = 1$, $\Lambda^L = 0$, $\Lambda^S = 0$, 
$H(T,P) = H_o^{LV}(T_V) + \int_{T_V}^T {C_p^{V}(\tau,P_{\rm ref})} d\tau
+ \int_{P_{\rm ref}}^P { [-T \alpha^V (T, p) +1] }dp$ ,\\ 
\end{center} 
with $\alpha^V$ the thermal expansion coefficient of vapor.

\vspace{2mm}
%\noindent
Since the heat capacity is strictly positive and 
$\Delta H^{fus}$ and $\Delta H^{vap}$ are also positive
(by our choice of $h_{\rm ref}$),
the dependence of $H$ on $T$ is monotonic, as shown schematically 
in Fig. \ref{fig:3}. Thus $T$ can be found from $H$. 

\begin{figure}[ht]   
 \centering
 \includegraphics[width=0.95\textwidth]{fig3} % ,height=3.0in
 \label{fig:3} \caption{
Schematic dependence of temperature $T$ and phase fractions 
on volumetric enthalpy $H$.}
\end{figure}


The PDE \eqref{eq:3:1} 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 $t_n$, 
we evaluate the thermophysical property values, compute the flux
$\vec{Q}$, and update the enthalpy $H$ from the discretized PDE 
to time $t_{n+1}$. 
From $H$ we find the phase indicators
${\bf{H_o^S}}$, ${\bf{H_o^L}}$, ${\bf{H_o^{LV}}}$, ${\bf{H_o^V}}$,
which determine the temperature and phase fractions as described 
above.


\subsection{Knudsen Layer and Vaporization temperature}\label{S:3:2}
%
Whereas melting occurs at a known melt temperature $T_m$ 
(= 1358 K for copper),
vaporization can occur at any temperature $T_V$, up to the critical
temperature $T_{\rm crit}$ ($\sim 8000$ K for copper), 
and must be determined at each time step.
$T_V$ 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. 
\cite{anisimov68,ytrehus77,knight79,knight81,miotelo95,ohmura98,
ellegaard99,gusarov05,gusarov00}.
A useful treatment is given in Gusarov-Smurov \cite{gusarov05}.
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 $T_{V}$, $P_{V}$ be the values at the melt surface 
(vaporization front),
and $T_{KL}$, $P_{KL}$ denote the values at the outer side of the
Knudsen layer (away from the target, towards the plume).
If the pressure ratio $P_{ratio} := P_{KL}/P_{V} < 1$ then 
subsonic evaporation occurs.  
If $P_{ratio} > 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 $T_{V}$,  
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 \cite{gusarov05}
are applied in an iterative procedure that determines $T_{V}$ 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 ${\bf H_o^V (T_V, P_{\rm ref})} < H < {\bf H_o^V (T_V, P_{\rm ref})}$)
at the current $T_V$. 
From the temperature $T_{KL}$, pressure $P_{KL}$, and gas density
$\rho^{KL}$ on the plume side, we find the
recession velocity $v_{V}$ of the vaporization front 
(by the Hertz-Knudsen relation and conservation of mass). 
From the recession velocity $v_{V}$ we update the vapor fraction of 
the (control) volume $Vol$ during a time-step $\Delta t$  via
$\Lambda^V$ = $\Lambda_{old}^V + v_{V} Vol / \Delta t$.
Then we estimate the new enthalpy of the volume as 
$H_{new}(T_V) = H_o^{LV}(T_V) + \Lambda^V \Delta H^{vap}(T_V)$,
which should equal $H$.  This constitutes an equation
$H_{new}(T_V) = H$ for the unknown $T_V$, 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 \cite{TEOS} for copper).


\section{Numerical Simulations}\label{S:4}

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 \S\ref{S:3}.

The target, of thickness 12 ${\mu}$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 $\Delta z = 5 \times 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 $\Delta t = 1 \times 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. \ref{fig:2}.
Thermodynamic properties of Cu are obtained from the ITTEOS 
\cite{TEOS} 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. \ref{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 ($T_b=2835K$) 
and vaporization becomes sigificant.
Thus the vaporization (surface) temperature $T_V$ and pressure $P_{V}$, 
as well as the pressure at the outer side of the KL, $P_{KL}$, 
increase quickly as seen in Fig. \ref{fig:6}--\ref{fig:7}.
Since $P_{KL} < P_{V}$ the evaporation is subsonic (see \S\ref{S:3:2}).

\begin{figure}[hpt]
 \centering
 \includegraphics[width=0.95\textwidth]{fig4}
% , height=4.0in]{fig4}
\label{fig:4}
 \caption{ Melt and vaporization depths vs time.}
\end{figure}

When the surface reaches temperatures near the critical state, 
the melt becomes transparent.  This bleaching of liquid properties 
occurs at temperatures about 0.9\,$T_{\rm crit}$ 
	\cite{batanov73}, \cite{karapetyan75}, \cite{yoo2000}.
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 $T_V \approx 0.9\,T_{\rm crit}$ as can be seen in 
Fig. \ref{fig:5}.

\begin{figure}[hpt]      
 \centering
 \includegraphics[width=0.95\textwidth]{fig5} %,height=3.8in
\label{fig:5}  \caption{Vaporization temperature vs time. }
\end{figure}


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
	\cite{zelldovich66}.
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 $\approx 16$ ns, 
the plasma is so dense that it completely shields the target. 
As a result, the vaporization temperature and pressure decrease,
and when $P_{V} < P_{KL}$ we arrive at the condensation regime 
(Fig. \ref{fig:6}).  During this stage, vapor from the plume 
recondenses on the target, the recession velocity becomes negative 
(Fig. \ref{fig:7}), and the crater depth decreases \eqref{fig:8}. 

\begin{figure}[hpt]   
 \centering
 \includegraphics[width=0.95\textwidth]{fig6} % ,height=4.0in
\label{fig:6}  \caption{ 
	Vaporization (saturation) pressure $P_{V}$ vs time, and 
	pressure at the outer side of the Knudsen layer, $P_{KL}$, 
	for comparison.
	}
\end{figure}

\begin{figure}[ht]   
 \centering
 \includegraphics[width=0.75\textwidth]{fig7} %,height=3.2in
\label{fig:7}  \caption{ 
	Vaporization (recession) velocity vs time.
	}
\end{figure}

After a while, at $\approx 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. \ref{fig:2}).

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 $\approx 0$ and constant evaporation 
depth of $\approx 150$ nm, as seen in Fig. \ref{fig:7} -- \ref{fig:8}.

\begin{figure}[ht]     
 \centering
 \includegraphics[width=0.75\textwidth]{fig8} % ,height=3in
\label{fig:8}  \caption{ 
	Vaporization depth vs time. Maximum depth reaches 200 nm shortly after 15 ns (when laser intensity peaks).
	}
\vspace{-2mm}
\end{figure}


\section{Conclusion}\label{S:5}

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 \S\ref{S: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 \S\ref{S:4}, 
the expected physical behavior of the system is captured very well.


\subsection*{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. \ref{fig:1}).
Support of the second author by the 
Flemish Fund for Scientific Research (FWO Vlaanderen)
is gratefully acknowledged.


\begin{thebibliography}{00}

\bibitem{book} Vasilios Alexiades and Alan D. Solomon,
\emph{Mathematical Modeling of Melting and Freezing Processes},
Hemisphere Publ.Co., 1993.

\bibitem{anisimov68} S. I. Anismov,
"Vaporization of metal absorbing laser radiation",
\emph{Sov. Phys., JETP}, \textbf{27}: 182-183, 1968.

\bibitem{batanov73} V. A. Batanov, F. V. Bunkin, A. M. Prokhorov, 
and V. B. Fedorov;
``Evaporation of metals caused by intense optical radiation'',
\emph{Sov. Phys., JETP}, \textbf{36}: 311-322, 1973.

\bibitem{becker02} J. S. Becker,
``Applications of inductively coupled plasma mass spectrometry and laser 
ablation inductively coupled plasma mass spectrometry in materials science",
\emph{Spectrochim. Acta, Part B}, \textbf{57}: 1805--1820, 2002.

\bibitem{bleiner05} D. Bleiner, 
``Mathematical modelling of laser-induced particulate formation in direct solid microanalysis",
\emph{Spectrochim. Acta Part B}, \textbf{60}(1): 49-64, 2005.

\bibitem{bleiner06} D. Bleiner and A. Bogaerts,
``Computer simulations of laser ablation sample introduction for 
plasma-source elemental microanalysis",
\emph{J. Anal. At. Spectrom.}, \textbf{21}: 1161--1174, 2006.

\bibitem{bogaerts03} Annemie Bogaerts, Zhaoyang Chen, Renaat Gijbels, Akos Vertes,
``Laser ablation for analytical sampling: what can we learn from modeling?",
\emph{Spectrochimica Acta B}, \textbf{58}: 1867-1893, 2003.

\bibitem{broekaert02} J. A. C. Broekaert,
\emph{Analytical Atomic Spectrometry with flames and plasmas}, 
Wiley VCH, 2002.

\bibitem{chen05} Zhaoyang Chen, Annemie Bogaerts;
``Laser ablation of Cu and plume expansion into 1 atm ambient gas",
\emph{J. Appl. Phys}, \textbf{97}: 063305, 2005.

\bibitem{chichkov96} B. N. Chichkov, C. Momma, S. Nolte, F. von Alvensleben, 
A. T\"{u}nnermann,
``Femtosecond, picosecond and nanosecond laser ablation of solids",
\emph{Appl. Phys. A}, \textbf{63}: 109--115, 1996.

\bibitem{chrisey94} D. B. Chrisey and G. K. Hubler, 
\emph{Pulsed Laser Deposition of Thin Films},
Wiley, 1994.

\bibitem{ellegaard99}
O. Ellegaard, J. Schou, H. M. Urbassek,
``Monte Carlo description of gas flow from laser-evaporated silver",
\emph{Applied Physics A}, \textbf{69}: S577­S581, 1999.

\bibitem{garrelie98} F. Garrelie, J. Aubreton, and A. Catherinot,
``Monte Carlo simulation of the laser-induced plasma plume 
expansion under vacuum: Comparison with experiments",
\emph{J. Appl. Phys.}, \textbf{83}, 5075, 1998.

\bibitem{gunther02} D. G\"{u}nther,
``Laser ablation-inductively coupled plasma mass spectrometry trends",
\emph{Anal. Bioanal. Chem.}, \textbf{372}: 31--32, 2002.

\bibitem{gusarov00} A. V. Gusarov, A. G. Gnedovets, and I. Smurov,
``Two-dimensional gas-dynamic model of laser ablation in an ambient gas,"
\emph{Applied Surface Science}, \textbf{154-155}: 66-­72, 2000.

\bibitem{gusarov05} A. V. Gusarov, I. Smurov,
``Thermal model of nanosecond pulsed laser ablation: Analysis of energy and mass transfer",
\emph{J. Applied Physics}, \textbf{97}: 014307, 2005.

\bibitem{harrington67} R. E. Harrington,
``Application of the Theory of Heat Conduction to the Absorption of 
Blackbody Radiation",
\emph{J. Applied Physics}, \textbf{38}: 3266, 1967.

\bibitem{itina97}
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",
\emph{J. Chem. Phys.}, \textbf{106}, 8905, 1997.

\bibitem{itina02} T. E. Itina, J. Hermann, Ph. Delaporte, and M. Sentis,
``Laser-generated plasma plume expansion: Combined continuous-microscopic 
modeling", \emph{Phys. Rev. E}, \textbf{66}, 066406, 2002.

\bibitem{itina03} T. E. Itina, J. Hermann, Ph. Delaporte, and M. Sentis,
``Combined continuous--microscopic modeling of laser plume expansion",
\emph{Appl. Surf. Sci.}, \textbf{27}: 208--209, 2003.

\bibitem{karapetyan75} R. V. Karapetyan and A. A. Samokhin,
``Influence of an increase in the transparency on the intense evaporation 
of metals by optical radiation",
\emph{Sov. J. Quantum Electron.}, \textbf{4}: 1141--1142, 1975.

\bibitem{knight79} C. J. Knight,
``Theoretical modeling of rapid surface vaporization with back pressure",
\emph{ AIAA Journal}, \textbf{17}(5): 519--523, 1979.

\bibitem{knight81} C. J. Knight,
``Transient Vaporization from a Surface into Vacuum",
\emph{AIAA Journal}, \textbf{20}(7): 950-954, 1981.

\bibitem{kools93} J. C. S. Kools,
``Monte Carlo simulations of the transport of laser-ablated atoms in a diluted gas",
\emph{J. Appl. Phys.}, \textbf{74}, 6401, 1993.

\bibitem{TEOS} P. R. Levashov and K. V. Khishchenko,
ITTEOS 5.8 software for calculation of EOS for metals, 2007.

\bibitem{lindner09} 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",
\emph{Anal. Chem.}, \textbf{81}: 4241--4248, 2009.

\bibitem{miotelo95} A. Miotello, A. Peterlongo, R. Kelly,
``Laser-pulse sputtering of aluminium: gas-dynamic effects with 
recondensation and reflection conditions at the Knudsen layer,"
\emph{Nuclear Instruments and Methods in Physics Research B}, 
\textbf{101}: 148-155, 1995.

\bibitem{noorbatcha87}
I. NoorBatcha, R. R. Lucchese, and Y. Zeiri,
``Monte Carlo simulations of gas-phase collisions in rapid desorption 
of molecules from surfaces",
\emph{J. Chem. Phys.}, \textbf{86}, 5816, 1987.

\bibitem{ohmura98} E. Ohmura, I. Fukumoto, I. Miyamoto,
``Molecular dynamics simulation of laser ablation of metal and silicon", 
\emph{Int. J. Jpn. Soc. Prec. Eng.}, \textbf{32}: 248-­253, 1998.
  Int. J. Jpn. Soc. Prec. Eng. 32

\bibitem{ready71} J. F. Ready,
\emph{Effects of High Power Laser Radiation}, 
Academic Press, 1971.

\bibitem{vertes93} A. Vertes, R. Gijbels, F. Adams, 
\textbf{Laser Ionization Mass Analysis}, 
Wiley VCH, 1993.

\bibitem{vonAllmen87} M. von Allmen,
\textbf{Laser-Beam Interactions with Materials},
Springer, 1987.

\bibitem{yoo2000} 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",
\emph{Appl. Phys. Lett.}, \textbf{76}: 783-785, 2000.

\bibitem{ytrehus77} T. Ytrehus,
``Theory and experiments on gas kinetics in evaporation",               
pp.1197-1212 in 
\emph{Rarefied Gas Dynamics},  edited by Potter, AIAA, 1977.

\bibitem{zelldovich66}
B. Zelldovich and Yu P. Raizer,
\emph{Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena },
Academic Press, 1966.

\end{thebibliography}


\end{document}