Numerical Simulation of Thermal Processes in a Domain of Thin Metal Film Subjected to an Ultrashort Laser Pulse

A thin metal film subjected to an ultrashort laser pulse is considered. With a sufficiently high laser intensity the process of the film heating may cause metal melting and even ablation. In this work, the numerical model of the melting and resolidification processes is presented. The mathematical model is based on the dual phase lag equation in which two positive constants appear, this means the relaxation and thermalization times. The considered equation contains a second-order time derivative and higher order mixed derivative in both time and space and should be supplemented by the appropriate boundary and initial conditions. The model of the melting and resolidification is presented in two versions. The first can be called ‘the introduction of the artificial mushy zone sub-domain’, while the second ‘the two forms of the basic energy equation’. At the stage of numerical computations, the implicit scheme of the finite difference method is used. The numerical algorithm is tested for the two proposed models which are applied to the computations concerning the thermal processes occurring in the cylindrical micro-domain (chromium, gold) subjected to an ultrashort laser pulse.


Introduction
In this paper, the application of the dual-phase lag equation (DPLE-see: Section 2) [1][2][3][4] for numerical modeling of the problems related to the microscale heat transfer are discussed. In particular, a thin metal film subjected to a laser pulse is analyzed. Generally speaking, the differences between the macroscopic heat conduction equation basing on the Fourier law and the models describing the ultrafast laser pulse interactions with the metal films appear because of an extremely short duration, extreme temperature gradients and very small dimensions of considered domain [5]. The delay times which are irrelevant in the case of the macroscopic heat transfer in the metals, play an important role in the analyzed case. As a rule, the DPLE of the first-order is considered (this is determined by the number of components in the Taylor series development), but the other forms of DPLE are also discussed. For example, both the heat flux q and the temperature gradient ∇T in the generalized Fourier law are expanded using the second-order Taylor formula (a second-order DPLE) [6,7]. The numerical solution of the second-order DPLE using the implicit scheme of the finite difference method (FDM) is reported by Chiriţă [8]. The mixed variants of DPLE are also considered, e.g., the second-order Taylor expression of heat flux and the first-order Taylor expression of temperature gradient [9]. The physical correctness of the higher order DPLE solutions and the limitations concerning the lag times are discussed by Fabrizio and Quintanilla [10,11]. This paper consists of the following parts. At first, the dual phase lag equation is presented and the analyzed problem is formulated. Next, the DPLE solution based on the implicit scheme of the finite difference method is described and the results of computations are shown. In the final part of the paper, conclusions are formulated.

Governing Equations
The start point for the further considerations is a generalized form of the Fourier law [1][2][3][4] q X, t + τ q = −λ∇T(X, t + τ T ), (1) where q is a heat flux vector, ∇T is a temperature gradient, λ is a thermal conductivity, X, t denotes the geometrical co-ordinates and time. The positive constants τ q , τ T correspond to the relaxation time and thermalization time, respectively. The relaxation time τ q takes into account the small-scale response in time, while the thermalization time τ T takes into account the small-scale response in space [43].
From Equation (2) results − q(X, t) = τ q ∂q(X, t) ∂t Now, the well-known thermal diffusion equation c ∂T(X, t) ∂t = −∇ · q(X, t) + Q(X, t) is used. In this equation c is the volumetric specific heat and Q (X, t) is the capacity of internal heat sources. Introducing Equation (3) to Equation (4) one obtains Because (c.f. Equation (4)) thus, in the end [1,2] c ∂T(X,t) ∂t The capacity of internal heat source Q(X, t) = Q(r, z, t) (the axisymmetric problem is analyzed, as shown in Figure 1) in the considered case is a sum of two components. The first component Q L (r, z, t) results from the laser action and according to the literature (e.g., [44]) the equation determining this value is the following where I 0 [J/m 2 ] is a laser intensity, t p [s] is a characteristic time of laser pulse, δ [m] is an optical penetration depth, R is a reflectivity of irradiated surface, r D [m] is a laser beam radius. The derivative of Q L with respect to time can be found analytically. The second component Q ph (r, z, t) complementary the source function Q(r, z, t) is associated with the phase transformations and will be described in the next section.

One-Domain Approach
The internal heat source resulting from the phase change is proportional to the local and temporary melting/solidification rate, in particular [45,46] where L is a volumetric latent heat, f S is a volumetric solid-state fraction in the neighborhood of the point considered, f L = 1 − f S . The function f L is a temperature dependent and for the border temperatures limiting the mushy zone sub-domain takes the values f L (T 1 ) = 0 and f L (T 2 ) = 1 (T 1 and T 2 correspond to the beginning and the end of the melting process). Thus where L is a volumetric latent heat. Finally, the source term controlling the evolution of latent heat can be written in the form Very often the course of the function f S between the border temperatures is assumed in the linear form, e.g., [46] One can see that this function fulfils the conditions f L (T 1 ) = 0 and f L (T 2 ) = 1. The first derivative of this function is equal to 1/(T 2 − T 1 ), while the second derivative is equal to 0 and the source function (11) takes a form Thus, the Equation (7) can be written as follows while the parameter C is called 'a substitute thermal capacity' and where c L , c S are the volumetric specific heats of liquid and solid state, while (for example) c M is the arithmetic mean of c L and c S . the dual phase la. Additionally, the following definition of the thermal conductivity is introduced where λ L , λ S are the thermal conductivities of liquid and solid state, while (for example) λ M is the arithmetic mean of λ L and λ S . It should be pointed out that Equation (14) containing the parameters C(T) and λ(T) defined according to (15) and (16) describes the thermal processes proceeding in a whole, conventionally homogeneous metal domain. In the case of the macroscopic melting/solidification models such an approach can be called 'a one-domain method'.
The mushy zone sub-domain limited by the border temperatures T 1 and T 2 appears in the case of the typical alloys melting/solidification. The melting or solidification of the pure metals proceeds at the constant temperature T*. To apply the one-domain approach for this case, the artificial mushy zone should be introduced. The artificial mushy zone is generated by 'a stretching' of solidification point T* to a certain interval [T* − ∆T, T* + ∆T] and for this interval the function f L is assumed in the linear form (as previously), this means and then The testing computations show that the assumed value of the interval ∆T (within a reasonable range of the few degrees e.g., 2.5-4.5 K) does not change significantly the results of numerical simulations [47].

Two Forms of the Dual Phase Lag Equation
The other approach to the numerical modeling of the melting/resolidification process discussed in this paper consists of the local application of two mutually exclusive forms of Equation (7). For the sub-domain in which the temperature did not reach the melting point, the dual phase lag equation with the source function (9) resulting from laser action is used, while the source function related to the melting and resolidification is equal to zero. The same situation takes place at the stage of the cooling process. On the other hand, in the sub-domain where the temperature is equal to the melting point, the source function related to the melting/resolidification process must be taken into account, while the temperature derivatives with respect to time are equal to zero.
Let us rewrite the Equation (7) with more precisely defined source function (9) controlling the melting/resolidification process c ∂T(r,z,t) ∂t Parameters c and λ differ for solid and liquid phases, of course. At the stage of cooling/heating process the function f L takes values 0 or 1 (solid or liquid) and the source function related to the melting/resolidification is equal to zero. Thus, for the subdomains in which the cooling/heating process takes place one obtains c ∂T(r,z,t) ∂t In turn, for the subdomains being at the stage of melting/resolidification process, the temperature is a constant value (corresponding to the solidification point T*) and then the cooling/heating rate is equal to zero, while the Equation (19) is of the form At the stage of numerical computations the implicit scheme of FDM is used. For each grid node the control volume (here in the form of rings) is assigned. Depending on the current situation the thermal processes in the volume considered are described by Equation (20) or (21).

Boundary and Initial Conditions
The cylindrical domain is limited by the side surface and two bases. The surfaces r = R 0 (external radius of domain) and z = Z (bottom base of cylinder) are far enough away from the source resulting from the laser's action, that adiabatic conditions can be accepted both for R 0 and Z. A similar condition can be assumed for z = 0. The laser action is taken into account by the introduction of an internal heat source, while the time of heat released outside (at the second stage of the process) is very short and this effect can be omitted.
The no-flux condition in the case of the first-order DPLE application is of the form [14,30] − λ n · ∇T(r, where n is the normal outward vector, n · ∇T(r, z, t) is the derivative of temperature in the normal direction. In the considered case the directional derivative corresponds to ∂T/∂r or ± ∂T/∂z. The initial temperature of the metal domain and the initial heating rate are also known where T p is the constant initial temperature.

Numerical Algorithm
To solve the problem discussed, the implicit scheme of the FDM is used. At first, the uniform time grid with the constant time step ∆t is introduced The axially symmetrical metal domain with dimensions R 0 and Z is covered by the regular geometrical mesh with step h, as shown in Figure 1.
The FDM algorithm will be described for the heating/cooling processes. The FDM equations for the melting/resolidification stages are very similar. The axially symmetrical metal domain with dimensions R0 and Z is covered by the regular geometrical mesh with step h, as shown in Figure 1. The temperature The FDM algorithm will be described for the heating/cooling processes. The FDM equations for the melting/resolidification stages are very similar. The following finite difference approximation of Equation (20) where (c.f. [32,45]) In the case of the axially symmetrical domain The following finite difference approximation of Equation (20) is proposed where (c.f. [32,45]) In the case of the axially symmetrical domain and The functions Φ can be called the shape functions of the FDM mesh, while R are the thermal resistances between the neighboring nodes, time level s = f or s = f − 1. The detailed mathematical considerations concerning the FDM approximation of the operator ∇(λ∇T) for different co-ordinate systems can be found in [45]. Thus, the Equation (25) can be written in the form Denoting one obtains The approximation of the no-flux boundary conditions (22) is assumed in the form Because the normal derivative for the bottom and upper surface is equal to ± ∂T/∂z, while for the cylinder side ∂T/∂r, consequently the approximation of n · ∇T(r, z, t) is very simple (the left and right sides differential quotients are used). Thus, the final form of conditions (33) will not be presented here. Obtained in this way the system of linear equations for transition t f − 1 → t f is solved using the iterative method. It should be noted that the presented algorithm is unconditionally stable [34].

Results of Computations
The calculations presented below are performed for a thin metal film of dimensions Z = 100 nm, R 0 = 100 nm made of chromium and subjected to a laser pulse (laser beam radius r D = R 0 /8, reflectivity R = 0.93, optical penetration depth δ = 15.3 nm, c.f. Equation (8)). The mathematical form of internal heat source resulting from the laser action suggests the orientation of the domain in the cylindrical coordinates (axially-symmetrical task). The dimensions Z and R 0 of the cylinder conventionally cut from the metal domain are assumed in such a way that the adiabatic conditions on the side and bottom surfaces can be accepted. The following values of chromium thermophysical parameters are assumed: thermal conductivities λ S = 93 W/(m K), λ L = 35 W/(m K), volumetric specific heats c S = 3.2148 MJ/(m 3 K), c L = 2.79276 MJ/(m 3 K) melting point T* = 2180 K, volumetric heat of fusion L = 2904 MJ/m 3 [48], relaxation time τ q = 0.136 ps, thermalization time τ T = 7.86 ps [1]. The solution obtained for the artificial mushy zone model corresponds to ∆T = 3 K, while the thermal conductivity and volumetric specific heat of this sub-domain are equal to the arithmetic means of the liquid and solid parameters. The initial temperature of the domain is equal to T p = 300 K. The step of the geometrical regular mesh h = 2 nm, while the time step ∆t = 0.0002 ps (the testing computations show that the further compaction of the grid does not change the numerical results).
The main purpose of the computations is to compare the melting/resolidification modeling results using two models, namely: model 1-melting at the temperature interval (artificial mushy zone), model 2-melting at the constant temperature.
The laser intensity and the characteristic time of the laser pulse are selected in such a way that the melting and resolidification take place, but the evaporation process does not appear. At the beginning, the laser intensity I 0 = 1.6 × 10 4 J/m 2 and three values of the characteristic times of laser pulse t p = 5 ps, t p = 10 ps, t p = 11 ps are taken into account (for t p > 11 ps the phase change has no place).
In Figure 2, the heating/cooling curves at points P 1 (0, 0) and P 2 (0, 20 nm) for t p = 10 ps are shown. At the most heated point P 1 the differences are visible, but at the other points (nodes) the courses of heating and cooling curves are almost identical. The next Figure shows the changes of the volumetric liquid state fraction at the point P 1 . The left part of Figure 3 illustrates the melting process, while the right part of this figure illustrates the resolidification process. Here, the differences are more visible, although the durations of these processes are comparable. The linear course of f L corresponds to the artificial mushy zone model (it results from theoretical considerations-Equation (12)).
Materials 2018, 11, x FOR PEER REVIEW 10 of 16 pulse tp = 5 ps, tp = 10 ps, tp = 11 ps are taken into account (for tp > 11 ps the phase change has no place). In Figure 2, the heating/cooling curves at points P1 (0, 0) and P2 (0, 20 nm) for tp = 10 ps are shown. At the most heated point P1 the differences are visible, but at the other points (nodes) the courses of heating and cooling curves are almost identical. The next Figure shows the changes of the volumetric liquid state fraction at the point P1. The left part of Figure 3 illustrates the melting process, while the right part of this figure illustrates the resolidification process. Here, the differences are more visible, although the durations of these processes are comparable. The linear course of fL corresponds to the artificial mushy zone model (it results from theoretical considerations-Equation (12)).   pulse tp = 5 ps, tp = 10 ps, tp = 11 ps are taken into account (for tp > 11 ps the phase change has no place). In Figure 2, the heating/cooling curves at points P1 (0, 0) and P2 (0, 20 nm) for tp = 10 ps are shown. At the most heated point P1 the differences are visible, but at the other points (nodes) the courses of heating and cooling curves are almost identical. The next Figure shows the changes of the volumetric liquid state fraction at the point P1. The left part of Figure 3 illustrates the melting process, while the right part of this figure illustrates the resolidification process. Here, the differences are more visible, although the durations of these processes are comparable. The linear course of fL corresponds to the artificial mushy zone model (it results from theoretical considerations-Equation (12)).   In Figure 4, the temperature histories at point P 1 (0, 0) for different values of characteristic times of laser pulse (for both models) are shown. Analysis of the mathematical form of the source function (8) confirms that for a smaller value of t p , the value of Q L (r, z, t) is greater, which causes that a higher local maximum of temperature is achieved. One can see that the time after which the local temperature reaches the maximum is significantly shorter.
The comparison of the obtained solutions shows, that the maximum differences between the temperatures determined by two considered models are equal to 325.5 K for t p = 5 ps, 51.5 K for t p = 10 ps and 17.1 K for t p = 11 ps. Smaller differences are observed for higher values of t p .
Materials 2018, 11, x FOR PEER REVIEW 11 of 16 In Figure 4, the temperature histories at point P1 (0, 0) for different values of characteristic times of laser pulse (for both models) are shown. Analysis of the mathematical form of the source function (8) confirms that for a smaller value of tp, the value of QL (r, z, t) is greater, which causes that a higher local maximum of temperature is achieved. One can see that the time after which the local temperature reaches the maximum is significantly shorter.
The comparison of the obtained solutions shows, that the maximum differences between the temperatures determined by two considered models are equal to 325.5 K for tp = 5 ps, 51.5 K for tp = 10 ps and 17.1 K for tp = 11 ps. Smaller differences are observed for higher values of tp. Similar computations are performed for the greater value of laser intensity, namely I0 = 2.6 × 10 4 J/m 2 and characteristic times of laser pulse tp = 19 ps, tp = 23 ps and tp = 28 ps, respectively (for tp < 19 ps the evaporation temperature is reached, while for tp > 28 ps the phase change has no place). In Figure 5 the temperature histories at the point P1 (0, 0) for different values of characteristic time of laser pulse and for both models are shown. The following differences are obtained: 149.2 K for tp = 19 ps, 75.7 K for tp = 23 ps and 8.3 K for tp = 28 ps. Similar computations are performed for the greater value of laser intensity, namely I 0 = 2.6 × 10 4 J/m 2 and characteristic times of laser pulse t p = 19 ps, t p = 23 ps and t p = 28 ps, respectively (for t p < 19 ps the evaporation temperature is reached, while for t p > 28 ps the phase change has no place). In Figure 5 the temperature histories at the point P 1 (0, 0) for different values of characteristic time of laser pulse and for both models are shown. The following differences are obtained: 149.2 K for t p = 19 ps, 75.7 K for t p = 23 ps and 8.3 K for t p = 28 ps. Comparing the two solutions presented above (Figures 4 and 5), one can notice that for the greater laser intensities and the greater characteristic times of laser pulse, the solution resulting from the artificial mushy zone model is closer to the results corresponding to the two forms of DPLE application.
The results of numerical computations allow, among others, to observe the temporary shape and dimensions of the pool of molten metal. In Figure 6 the final shape of the pool before the stage of resolidification is presented.  Comparing the two solutions presented above (Figures 4 and 5), one can notice that for the greater laser intensities and the greater characteristic times of laser pulse, the solution resulting from the artificial mushy zone model is closer to the results corresponding to the two forms of DPLE application.
The results of numerical computations allow, among others, to observe the temporary shape and dimensions of the pool of molten metal. In Figure 6 the final shape of the pool before the stage of resolidification is presented. Comparing the two solutions presented above (Figures 4 and 5), one can notice that for the greater laser intensities and the greater characteristic times of laser pulse, the solution resulting from the artificial mushy zone model is closer to the results corresponding to the two forms of DPLE application.
The results of numerical computations allow, among others, to observe the temporary shape and dimensions of the pool of molten metal. In Figure 6 the final shape of the pool before the stage of resolidification is presented.  To investigate more accurately the 'action' of the model using two forms of DPLE, the material with significantly different (in comparison with chromium) thermophysical parameters is taken into account. Thus, the cylindrical domain (R 0 = 100 nm, Z = 100 nm) made of gold is considered. The material parameters are the following: λ S = 315 W/(m K), λ L = 105 W/(m K), volumetric specific heats c S = 2.4897 MJ/(m 3 K), c L = 2.8757 MJ/(m 3 K) melting temperature T* =1336 K, volumetric heat of fusion L = 1229.99 MJ/m 3 . Relaxation time τ q = 8.5 ps, thermalization time τ T = 90 ps [49]. The laser beam parameters: r D = R 0 /8, laser intensity: I 0 = 5 × 10 4 J/m 2 , characteristic time of laser pulse: t p = 5 ps. The initial temperature of domain is equal to T p = 300 K.
In Figure 7, the heating/cooling curves at the selected points from the metal domain are shown.
The results of computations are (from a qualitative point of view) similar to those obtained previously. It can be seen that the maximum temperature at the node (0, 0) is obtained after the time of 11 ps, while in the case of chromium, this time is equal to 13.5 ps (for the same characteristic time of laser pulse). This may be surprising because chromium should be heated to a much higher temperature to initiate the melting process, which should last longer. It results from the fact that the thermal conductivity of the gold is significantly bigger, which causes the some equalization of the temperature field in the domain and significantly reduces local temperature gradients (heat fluxes). Additionally, the lag times of gold are much longer than in the case of chromium. These factors mean that despite such different materials, the heating/cooling curves with respect to time are quite similar. The local and temporary temperatures are, of course, quite different. To investigate more accurately the 'action' of the model using two forms of DPLE, the material with significantly different (in comparison with chromium) thermophysical parameters is taken into account. Thus, the cylindrical domain (R0 = 100 nm, Z = 100 nm) made of gold is considered. The material parameters are the following: λS = 315 W/(m K), λL = 105 W/(m K), volumetric specific heats = 2.4897 MJ/(m 3 K), = 2.8757 MJ/(m 3 K) melting temperature T* =1336 K, volumetric heat of fusion L = 1229.99 MJ/m 3 . Relaxation time τ = 8.5 ps, thermalization time τ = 90 ps [49]. The laser beam parameters: = R0/8, laser intensity: I0 = 5 × 10 4 J/m 2 , characteristic time of laser pulse: tp = 5 ps. The initial temperature of domain is equal to Tp = 300 K.
In Figure 7, the heating/cooling curves at the selected points from the metal domain are shown.
The results of computations are (from a qualitative point of view) similar to those obtained previously. It can be seen that the maximum temperature at the node (0, 0) is obtained after the time of 11 ps, while in the case of chromium, this time is equal to 13.5 ps (for the same characteristic time of laser pulse). This may be surprising because chromium should be heated to a much higher temperature to initiate the melting process, which should last longer. It results from the fact that the thermal conductivity of the gold is significantly bigger, which causes the some equalization of the temperature field in the domain and significantly reduces local temperature gradients (heat fluxes). Additionally, the lag times of gold are much longer than in the case of chromium. These factors mean that despite such different materials, the heating/cooling curves with respect to time are quite similar. The local and temporary temperatures are, of course, quite different.

Conclusions
Two variants of numerical modeling of melting and resolidification occurring in thin metal films subjected to an ultrashort laser pulse have been discussed. In the first variant, the artificial mushy zone was introduced (melting and resolidification in the interval temperature), while the second variant was based on two forms of a dual phase lag equation. From the numerical point of view, the developed computer programs are of similar complexity, although the program for the artificial mushy zone approach is somewhat simpler.
The discussion of the obtained results has been presented in detail in Section 7. Summing up, in the case of microscale heat transfer modeling the numerical solutions based on the two proposed models to some extent differ from each other, but it seems that in the practical applications such differences are acceptable. The closer results of solutions have been obtained for the greater values of Figure 7. Temperature history at the points P 1 (0, 0), P 2 (0, 20 nm), P 3 (0, 40 nm), P 4 (0, 60 nm) and P 5 (0, 80 nm), gold, I 0 = 5 × 10 4 J/m 2 , t p = 5 ps.

Conclusions
Two variants of numerical modeling of melting and resolidification occurring in thin metal films subjected to an ultrashort laser pulse have been discussed. In the first variant, the artificial mushy zone was introduced (melting and resolidification in the interval temperature), while the second variant was based on two forms of a dual phase lag equation. From the numerical point of view, the developed computer programs are of similar complexity, although the program for the artificial mushy zone approach is somewhat simpler.
The discussion of the obtained results has been presented in detail in Section 7. Summing up, in the case of microscale heat transfer modeling the numerical solutions based on the two proposed models to some extent differ from each other, but it seems that in the practical applications such differences are acceptable. The closer results of solutions have been obtained for the greater values of characteristic times of laser pulse t p . The heating/cooling curves presented in the previous section concerned the most intense regions of the heat exchange, while in the peripheral sub-domains, the differences between the solutions were insignificant.
The authors' approach to phase change modeling was not used. Taking into account the microtechnologies based on the melting or solidification phenomena, one can believe that this type of numerical models should be useful at the stage of technology design.
Further research in this area is expectd to tackle the phenomenon of ablation caused by the heating of material to very high temperatures.