AN AXISYMMETRIC MODEL FOR SOLID-LIQUID-VAPOR PHASE CHANGE IN THIN METAL FILMS INDUCED BY AN ULTRASHORT LASER PULSE

An axisymmetric model for thermal transport in thin metal films irradiated by an ultrashort laser pulse was developed. The superheating phenomena including preheating, melting, vaporization and re-solidification were modeled and analyzed. Together with the energy balance, nucleation dynamics was employed iteratively to track the solid-liquid interface and the gas kinetics law was used iteratively to track the liquid-vapor interface. The numerical results showed that higher laser fluence and shorter pulse width lead to higher interfacial temperature, larger melting and ablation depths. A simplified 1-D model could overestimate temperature response and ablation depth due to the omission of radial heat conduction.


INTRODUCTION
The great application potential of ultrashort lasers has opened an interesting research area in the micro-and nanoscale regime due to their ability to deposit extremely high density energy into a target in a very short period of time. Many research works have demonstrated that Fourier's Law of heat conduction is not suited for modeling ultrafast heat transport processes, in which the energy excitation times are shorter than or comparable to the relaxation times of heat carriers (Wang and Prasad, 2000). For that reason, the two temperature models have been employed for describing ultrashort laser-metal interactions, including three consecutive stages: 1) deposition of laser energy to electrons in a subsurface layer, 2) transport of energy to deeper part of electrons, and 3) transfer of energy from electrons to lattice (Rethfeld et al., 2002;Hohlfeld et al., 2000). Anisimov et al. (1974) originally proposed the two temperature model which was rigorously derived later by Qiu and Tien (1993) based on Boltzmann transport equation. The revised model considered hyperbolic heat conduction in electrons but neglected heat conduction in the lattice. Chen and Beraun (2001) proposed a dual-hyperbolic model for materials in which heat conduction in the lattice was taken into account. For the case of constant material properties, those two-temperature models can be represented by the dual-phase-lag models (Tzou, 1997;Tzou, 2006).
Most of the existing models for ultrafast solid-liquid and liquidvapor phase changes dealt with either pure conduction (Diniz Neto and Lima, 1994;Chen et al., 2002) or one-dimensional (1-D) melting or melting and vapor phase change problems (Zhang and Chen, 2007;Huang et al., 2009;Zhang and Chen, 2008). Under high laser fluence and short pulse width, melting could take place followed by solidification as heat diffused away. The rapid phase change induced by ultrashort pulsed lasers is controlled by nucleation dynamics at the interface, not by the interfacial energy balance (Von Der Linde et al., 1987). Recently, the authors (Zhang and Chen, 2008) proposed an implicit, fixed grid interfacial tracking method to solve kinetics controlled rapid melting and resolidification in a free-standing metal film. For a laser heating at high fluence, the surface temperature may exceed the saturation temperature and thus, vaporization takes place. Shi et al. (2007) analytically investigated solid-liquid-vapor phase changes of metal particles caused by a nanosecond laser. By assuming a fixed velocity of liquid-vapor interface and neglecting the ablation depth, Chowdhury and Xu (2003) proposed a numerical model to simulate the melting-vaporization-resolidification process during femtosecond laser-metal interaction.
The 1-D models for solid-liquid-vapor phase change (Huang et al., 2009) are valid only if the laser spot size is much larger than the target thickness. If the above condition is not met, a multi-dimensional model is needed. The authors recently proposed an axisymmetric model for melting and resilification (Baheti et al., 2010). To accurately describe the entire thermal process of an ultrashort pulsed laser interaction with a metal film, vaporization should also be taken into account if existing.
This work presents an axisymmetric model of ultrashort pulsed laser-metal interaction incorporating preheating, melting, vaporization, removal of vaporized material, and re-solidification. The numerical model is formulated based on a finite volume method. The solid-liquid and liquid-vapor interfaces are tracked and studied under different laser parameters. The vapor generated at the surface is assumed to be immediately removed so that the computational domain includes only solid and liquid phases. Together with the energy balance, temperature and velocity, controlled by nucleation dynamics for the solid-liquid interface and by the gas kinetics law for the liquid-vapor interface, are iteratively solved through an iterative computational procedure. The temperature-dependent thermophysical properties are incorporated for more accurate prediction of thermal transport in electrons and lattice over the time history.

METHODOLOGY
The problem under consideration here consists of a cylindrical gold sample of radius M and thickness L impinged by a Gaussian laser beam

Frontiers in Heat and Mass Transfer
Available at www.ThermalFluidsCentral.org of characteristic radius of r 0 , much smaller than the radius M. The central line of the incoming laser beam is assumed to coincide with the central line (r = 0) of the sample. The heating effect caused by the laser beam is axisymmetric in nature. Therefore, numerical analysis is performed only in the r and z directions.
The two dimensional axisymmetric energy equations for electrons and lattice are: e  e  e  e  l   T  T  T Equation (1) is valid in the entire computational domain, while Eq. (2) is valid for both solid and liquid phases but not at the solid liquid interface.
The time t = 0 is defined as the time when the pulse energy reaches its peak and the numerical simulation starts from t = -2t p . Therefore, the volumetric heat source of a Gaussian laser is given by: where r 0 is the spot characteristic radius of a Gaussian beam defined at e -1 position, β = 4ln(2), and the factor  is to correct the film thickness effect. Electron-lattice coupling factor denoted by G in Eqs. (1) and (2) is temperature-dependent (Chen et al., 2003).
where G RT is the coupling factor at room temperature, and A e and B l are constants. The value of G is taken 20% higher in liquids than that in solids as the electrons are more likely to collide in a liquid phase (Kuo and Qiu, 1996). The temperature-dependent electron heat capacity is approximated as (Chen et al., 2006 where C΄ is given as: Although most of the conduction in pure metals (about 99% for gold) is due to free electrons and the rest of about one percent is contributed by lattice at equilibrium (Klemens and Williams, 1986), lattice heat conduction is included in this work since the lattice energy balance at the solid-liquid interface plays an important role [see Eq. (11)] on melting and solidification.. The thermal conductivity of the electrons widely employed in two temperature models is given by: which is valid only when electron temperature is much lower than Fermi temperature. A more general form of electron thermal conductivity was proposed by Anisimov and Rethfeld (1997) where ζe = T e /T F and ζl = T l /T F . The two constants in Eq. (8) are χ = 353 W/m-K and η = 0.16 for gold. For higher electron temperatures, ζe >> 1, Eq. (8) results in the well known dependence k e ≈ (T e ) 5/2 , that is the characteristics of low density plasma. On the other hand, Eq. (10) reduces to Eq. (7) under low temperature limit, ζe << 1. The initial temperature (T i ) of the metal sample is: The heat loss due to radiation from the surfaces is negligible as reported (Huang et al., 2009). All the boundaries of the sample are thus considered to be adiabatic: By expressing the position of solid-liquid interface as ( , ) z s r t s =  , energy balance at the interface is given by (Faghri and Zhang, 2006): where s u  is the solid-liquid interfacial velocity during melting or solidification.
The velocity of solid-liquid interface is governed by the energy balance specified by Eq. (11) if the melting process is conventional. For a rapid heating, although Eq, (11) still needs to be satisfied, the solidliquid phase change is dominated by nucleation dynamics and the solidliquid interface velocity is described as (Chen et al., 2006): The interfacial temperature ( ,I T  ) may be higher than the normal melting point during melting and lower during the solidification. For the liquid-vapor interface, energy balance at the interface together with the vaporization rates derived from the kinetics laws are applied to solve its temperature, velocity and location. Let the shape of the liquid-vapor interface be expressed as ( , ) v z s r t =  , the energy balance at an liquid-vapor interface may be given by (Faghri and Zhang, 2006): where v T  and v u  are the liquid-vapor interfacial temperature and velocity, respectively. The latent heat of vaporization h ν is related to T ℓν through (Xu et al., 1999): Specific heat (J/kg-K) 105.1 + 0.291T l -8.713 x 10 -4 T l 2 + 1.187 x 10 -6 T l 3 -7.051 x 10 -10 T l 4 + 1.538 x 10 -13 T l 5 (Chen and Beraun, 2001) 163.205 Thermal conductivity at equilibrium (W/m-K) 320.973 -0.0111T l -2.747 x 10 -5 T l 2 -4.048 x 10 -9 T l 3 (Klemens and Williams, 1986) 37.72 + 0.0711T l -1.721 x 10 -5 T l 2 + 1.064 x 10 -9 T l where h ν0 is the latent heat of vaporization at absolute zero, and T c is the critical temperature. With a known liquid-vapor interfacial velocity and h v given by (14), Eq. (13) can be used to solve the liquid-vapor interfacial temperature T ℓν .
To determine the liquid-vapor interfacial velocity, Clausius-Clapeyron equation is employed to describe the slope of saturation pressure-temperature curve assuming the process is under thermal equilibrium and ideal gas conditions: The pressure at the liquid-vapor interfacial temperature can be obtained by integrating Eq. (15): The molar evaporation flux (j ν ) at the surface can be calculated by the Hertz-Knudsen-Langmuir equation derived from the kinetic theory of gases (Birks et al., 2006): where A is an accommodation coefficient that shows the portion of vapor molecules striking the liquid-vapor surface is absorbed by this surface (Akhatov et al., 2001). Xu et al. (1999) recommended a value of 0.82 for this coefficient. With the pressure given by Eq. (16), the liquid vapor interfacial velocity v u  can be obtained from j ν by the relation:

NUMERICAL SOLUTION
In this work the five stages of thermal transport in a metal film with solid-liquid-vapor phase change caused by ultrafast laser heating include: 1. Pure heat conduction takes place in both electrons and lattice, described by Eqs. (1) and (2) with an internal laser heat source given by Eq. (5). 2. Liquid is formed from the heated surface as melting starts, leading to a solid-liquid interface moving toward the inner part of the film. The solid-liquid interfacial temperature and velocity are governed by Eqs. (11) and (12).  18). The liquid is considered to be removed immediately from the computational domain once it becomes vapor. In this stage both the solid-liquid and the liquidvapor interfaces are tracked simultaneously. 4. Evaporation ceases but melting still continues as the electrons and lattice reach equilibrium. The melting process gradually slows down after it reaches its maximum melting depth. 5. Liquid metal starts to cool down, leading to resolidification.
The computational domain is uniformly divided by 2D grids (Fig.  1). Each control volume is represented by a grid point at its center. During the liquid-vapor phase change, once the interface moves over a grid point, the respective grid point is considered under the vapor phase. This grid point (Fig. 2) is treated with a block-off method (Patankar, 1980), i.e. setting the k e and k l in the vapor zones to be zero. Accordingly, the control volume is changed depending on whether the vapor phase in a control volume is above or below the grid point (Fig.  2). When the liquid-vapor interface is below the grid point, the control volume is removed and the remaining liquid part is attached to the adjacent (right side) control volume as show in Fig. 2. Otherwise, the control volume is retained with the remaining liquid part.
The implicit finite volume equations are obtained by integrating Eqs. (1) and (2) in each control volume and time step (Patankar, 1980). Discretizing and simplification Eqs. (1) and (2) respectively over the control volume gives: where T ξ,P , T ξ,E , T ξ,W , T ξ,N and T ξ,S are respectively electron (ξ = e) and lattice (ξ = l) temperature of the grid points P, E, W, N and S (Fig. 1) where S p and G P are the heat source intensity and the electron-phonon coupling factor at grid point P; ,P o T ξ represents the electron (ξ = e) and lattice (ξ = l) temperature of grid point P at the previous time step; k ξ,e , k ξ,w , k ξ,n and k ξ,s are the thermal conductivities at the faces of control volume e, w, n and s ( Fig. 1) respectively and are obtained by harmonically averaging the thermal conductivities at the two adjacent control volumes (Patankar, 1980). The solid-liquid interfacial velocity can be obtained from the volumetric enthalpy and the energy equation for the control volume that contains the solid liquid interface (Zhang and Chen, 2007) where ,I o l T is the solid-liquid interfacial temperature at the previous time step. The third and the fourth terms in the bracket at the right hand side of Eq. (21) represent the effect of electron-lattice interaction and change of the interfacial temperature. The last term in the equation incorporates the effect of lateral direction on the velocity given by Eq. (11). The liquid fraction f P in the control volume P (considering axial direction) that contains interface is: Equation (21) is used along with Eq. (12) to determine the solid-liquid interfacial velocity and temperature. The solid-liquid interfacial location is then determined by using: Similarly, the liquid-vapor interfacial location is updated by: The numerical solution starts from time t = -2t p . Before the onset of melting, the electron and lattice temperatures are obtained by solving the coupled equations in Eq. (19) with ξ = e and l using an alternating direction implicit (ADI) method. The numerical solution is performed for the whole computational domain in a layer by layer pattern, incorporating the heat flow in both rand zdirections, i.e. the grids of the first column (r = 0) is solved from the top to the bottom, following the second and so on until the last column (r = M). Melting starts from the control volume with r = z = 0. Since the rapid melting and re-solidification are controlled by nucleation dynamics, the interfacial temperature is unknown and is related to the interfacial velocity by Eq. (12). For those control volumes with partial melting, an iterative procedure is employed to solve the interfacial temperature, velocity and location at each time step. For the time step when the lattice temperature of a control volume first exceeds the melting point, the lattice temperature in that control volume is first computed by assigning the initial values T l,P = T m , ,P l a =10 20 , and l b = ,I l T x10 20 (Zhang and Chen, 2007). The iterative procedure is as follows: 1. The solid-liquid interfacial temperature, T l,I , is assumed and Eq. (21) is used to determine the solid-liquid interfacial velocity. The liquid fraction and location in the first control volume are determined using Eqs. (22) and (23), respectively. 2. Nucleation dynamics, Eq. (12), gives a solid-liquid interfacial velocity which is then compared with the interfacial velocity obtained from Step 1. Depending on whether the interfacial velocity obtained from Eq. (21) is higher or lesser than that obtained from Eq. (12), the interfacial temperature will be increased or decreased by an appropriate value. 3. The two equations in Eq. (19) are solved simultaneously to obtain the electron and lattice temperatures in the whole domain. 4. Steps 1 to 3 are repeated until the maximum difference between the two consecutive interfacial velocities obtained from eq. (24) and (12) is less than 10 -3 m/s. 5. The interfacial location for the second layer is updated using Eq.
(22) with the interfacial velocity obtained from step 4. The process repeats until the last layer is solved.
Similarly, the iterative procedure is applied to solve the liquidvapor interfacial temperature, velocity and location for those control volumes with partial vaporization. For the time step when the lattice temperature of a control volume first exceeds the saturation temperature, the lattice temperature in that control volume is set at the saturation temperature. The iterative procedure is as follows: 1. Assume an interfacial velocity * (24) where α ℓν is the appropriate under-relaxation factor used to obtain a converged solution. In this paper, its value is set to be 0.8. 6. Update * v u  with the value of *** v u  . 7. Repeat Steps 2 to 6 until the difference between the interfacial velocities obtained from the consecutive iteration is less than a predefined precision value, 10 -5 m/s in this work.

RESULTS AND DISCUSSIONS
A gold film sample of 0.55-μm in radius and 0.55-μm in thickness was considered for numerical analysis. The initial temperature of the sample was taken to be 300 K. Convergence of model mesh and time step were first studied with three uniform grid meshes, consisting of 25×25, 50×50, 100×100 grids in the computational domain. It was found that the model of 50x50 grids with a time step of tp/200 yielded a reasonably convergent solution, and thus they were used in the following simulations. The validity of the present axisymmetric model was then verified with the 1-D model (Huang et al., 2009). To compare the present results with the 1-D model, a uniform laser beam was applied on the entire top surface of the axisymmetric model by setting r0 in Eq. (5) as infinity. The thermophysical properties of gold used for the comparison are listed in Table 1, which are the same as those used in our former work (Huang et al., 2009). Figure 3 compares the liquid-vapor interfacial temperature, velocity and location between the present axisymmetric 2-D model and the 1-D model for different laser fluences. The laser pulse width used in the calculations was t p = 100 fs. The liquid-vapor interfacial temperature reaches its peak value 4785 K at 23 ps for the pulse of Jo = 0.7 J/cm 2 and 4040 K at 23 ps for the pulse of 0.6 J/cm 2 , obtained by the present axisymmetirc model. They are almost identical to the 1-D results. The time histories of the interfacial velocities in Fig. 3(b) and ablation depth in Fig. 3(c) also show the excellent agreement between the two models. Further comparison was conducted for a longer laser pulse of t p = 2.5 ps. Again, the results agree with each other very well. Those results are not presented here for brevity. Figure 4 shows the surface lattice temperature distributions at t = 2 ps for fluences Jo = 0.65 J/cm 2 and 0.75 J/cm 2 . The simulations were performed for the laser pulse of tp = 100 fs and r0 = 0.35 μm. Figure 5 plots the lattice temperature contours at the time when the lattice temperatures reached their maximum, 3960 K at t = 20.3 ps for Jo = 0.65 J/ cm 2 and 4658 K at t = 21.2 ps for 0.75 J/ cm 2 . Interaction with a 500-fs pulse was also investigated. As expected, the peak lattice temperatures are slightly lower, 3896 K and 4566 K respectively. It is noted that at these two fluences evaporation hardly happens for laser pulses longer than 500 fs.
(a) t p = 100 fs, J 0 = 0.65 J/cm 2 (b) t p = 100 fs, J 0 = 0.75 J/cm 2 Fig. 4 Comparison of the top surface lattice temperature at t = 2 ps for different fluences Distributions of lattice temperature along the centerline (r = 0) at t = 2 ps and 10 ps are depicted in Fig. 6 for two laser pulses of 100 fs and 500 fs. Figure 7 shows the maximum ablation depths. For the 500-fs pulse, the maximum ablation depths were 0.0015 nm for 0.65 J/ cm 2 and 0.0062 nm for 0.75 J/ cm 2 , which are slightly smaller than those for the 100-fs pulse, 0.0017 nm and 0.0067 nm. Figures 8(a)-(c) tracks the solid-liquid interfacial temperature, velocity and location at the center of the film irradiated by the 100-fs pulse. As shown in Fig. 8(a), the maximum lattice temperature is below the saturation temperature for the laser fluences 0.50 J/ cm 2 and 0.30 J/ cm 2 . The normal boiling point of gold at 1 atm is 3127 K. No vaporization occur in these two cases. For a higher fluence of 0.75 J/ cm 2 , for example, the solid-liquid interfacial temperature reaches 4389 K for the case of 0.75 J/ cm 2 as shown in Fig. 8(a). It is noted that the interfacial temperature reaches their maximum value at approximately 20 ps, which is consistent with the results reported in literature (Zhang and Chen, 2007;Jing et al., 2009). Figure 8(b) shows the change in solid-liquid interfacial velocity versus time is similar to lattice temperature, rising in the early time (about 20 ps) and then gradually decreases. The solid-liquid interfacial locations, as functions of time are shown in Fig. 8(c). Apparently, the times when the peak melting depths appear differ from those when the maximum lattice temperatures are present. After reaching the peak melting depth, resolidification starts to take place. The metal film completely resolidifies at approximately 241 ps and 845 ps for the laser fluence 0.30 J/ cm 2 and 0.50 J/ cm 2 , respectively. Figure 8(d) compares the melting depths at r = 0 and r = 100 nm. It is shown that the melted zones decreases as the radial location goes farther from the centre of the target. At t = 300 ps, the melting depths at r = 100 nm is 82 nm for 0.65 J/ cm 2 and 102 nm for 0.75 J/ cm 2 , compared to 91 nm and 110 nm at r = 0 respectively. The differences of the melting depths are similar at t = 1 ns. Liquid-vapor interfacial temperature, velocity and location at different radius locations are respectively compared in Fig. 9 for the 100-fs laser pulse. Figure 9(a) shows that the interfacial temperatures at r = 100 nm also exceed the normal boiling point. As seen in Fig. 9(b), the trend of change in the liquid-vapor interfacial velocity is similar to that of the interfacial temperature. The maximum velocities occur at about t = 20 ps, same as the time when the solid-liquid interfacial velocities reach its maximum. Figure 9(c) illustrates the ablation depth. The maximum ablation depth at the centre (r = 0) is 0.0066 nm for 0.75 J/ cm 2 occurring at about 105 ps, while it is only 0.0015 nm at r = 143 nm. Comparing the results in Figs. 8(b) and 9(b) reveals that the evaporation rates are three order lower than the melting rates. It should be noted that ablation may not really occur at this fluence level since the simulated ablation depth is small compared to lattice constants. To demonstrate the effect of radial thermal diffusion on the superheating process, the volumetric laser heat source at the center of a Gaussian laser beam, i.e., r = 0 in Eq. (5) is applied for the 1-D model. Figures 10 and 11 compare the 1-D results with those at r = 0 of the present axisymmetric 2-D results. It can be seen clearly in Fig. 10 that the 1-D model predicts greater temperature and vaporization rate at the liquid-vapor interface than the axisymmetric model that includes the radial heat conduction effect. Figure 11(a) shows the comparison of maximum temperature for various fluences. Due to the fact that the 1-D model excludes radial heat conduction, it overestimates the temperature response by about 13%. The same tendency is shown for ablation depth; however, the differences can be over 100%.

CONCLUSIONS
An axisymmetric interfacial tracking method for solid-liquid and liquidvapor interfaces is developed to model melting, vaporization and resolidification in thin metal films irradiated by an ultrashort short Gaussian laser beam. Together with energy balance, nucleation dynamics is employed iteratively to track the solid-liquid interface and the gas kinetics law is used iteratively to track the liquid-vapor interface. Vaporized (ablated) material is removed from the computational domain by employing the block-off method once vaporization takes place. The computer code is first validated with the 1-D results. Numerical results show that higher laser fluence and shorter pulse width lead to deeper ablation depth, higher liquid-vapor interfacial temperature and greater interfacial velocity. The same trends are also found for the solid-liquid interface. It is also found that under the same laser heating conditions, the ablation depth is four orders smaller than the melting depth and the interfacial liquid-vapor velocity is three orders lower than solid-liquid interfacial velocity. Due to the fact that radial heat conduction is not taken into accounted, a simplified 1-D model could overestimate temperature response by about 13% and ablation depths by more than 100%.