Possible electronic entropy-driven mechanism for non-thermal ablation of metals

The physical mechanism for metal ablation induced by femtosecond laser irradiation was investigated. Results of calculations based on finite-temperature density functional theory (FTDFT) indicate that condensed copper becomes unstable at high electron temperatures due to an increase of electronic entropy at large volume, where the local density of states near the Fermi energy increases. Based on these results, an electronic entropy-driven (EED) model is proposed to explain metal ablation with a femtosecond laser. In addition, a mathematical model is developed for simulation of the laser ablation, where the effect of the electronic entropy is included. This mathematical model can quantitatively describe the experimental data in the low-laser-fuence region, where the electronic entropy effect is determined to be especially important.

Specific phenomena have been observed by irradiation of a metal surface with a femtosecond laser, such as ultrafast structural changes, 1,2 bond hardening, 3 and the emission of excessively high energy ions and neutral atoms. [4][5][6] The process of removing materials with an intense laser is called ablation, and ablation which cannot be explained under the assumption of thermal equilibrium such as the last one is referred to as non-thermal ablation. One of the reasons for the specific attention to non-thermal ablation from industry is that it can decrease the thermal damage region. [7][8][9] However, a complete understanding of non-thermal ablation of metals is still missing, so that optimal laser conditions for precision processing with femtosecond laser irradiation cannot be predicted from simulation.
Femtosecond laser irradiation on a metal surface changes the electron subsystem of the metal from the ground state into excited states. An electron subsystem is thermalized to the Fermi-Dirac distribution with electron temperature T e , via electron-electron interaction, of which the scattering time τ ee is approximately 10-100 fs in metals. 10,11 At the same time, the lattice temperature T l begins to increase by energy transfer from the electron subsystem via electron-phonon interaction, of which the scattering time τ el is on the time scale of picoseconds. [12][13][14][15] Therefore, under the assumption of instantaneous and local thermalization in the electron subsystem, T e ≫ T l is expected to be realized long before τ el by intense femtosecond laser irradiation. This is the main concept of the two-temperature model (TTM). 16 The TTM has been widely employed to simulate 2,3,17-20 such phenomena induced by femtosecond laser irradiation and it has been successful for quantitative description of the experimental data. 2,3 These TTM calculation results indicated that change of the interatomic forces due to high T e is required to reproduce the a) Electronic mail: tanaka@cms.phys.s.u-tokyo.ac.jp experimental data. Recently, some calculation studies have suggested that ablation can be caused by a change of the electronic states without emitting electrons from a metal surface. 18,19 This explanation, which does not require the neutrality breakdown, is supported experimentally, 21 although there is a conflict with a previous explanation based on the Coulomb explosion (CE) model, [4][5][6] which is verified in an insulator 22 and a molecule. 23 Therefore, in case of metals, where many electrons with good mobility exist, we expect that the explanation based on the CE model should be replaced by the former explanation. Although these previous studies have suggested that a large contribution of the force induced by change of the electronic states originates from the kinetic energy of free (delocalized) electrons, 18,19,24,25 the validity of this explanation has yet to be clarified. Furthermore, comparison with experimental results has not been performed sufficiently.
Here, we first investigate the physical mechanism of non-thermal ablation of copper (Cu) based on results of first-principles calculations, and propose the electronic entropy-driven (EED) model to describe it. Subsequently, a mathematical model is developed for simulation of the ablation depth, where the EED model is employed. Finally, we present results of the simulation and compare them with the experimental data. 26,27 It should be noted that the TTM is employed in all our calculations.
To investigate the main contribution that causes nonthermal ablation, we conducted first-principles calculations based on finite-temperature density functional theory (FTDFT). 28 Before thermal equilibrium is achieved, a system irradiated with an intense femtosecond laser can be represented as T e > T l by employing the TTM. To examine T e dependence of the stability of Cu, the electronic free energy F was calculated, which is defined as: where E is the internal energy and S is the electronic entropy. S for independent particles that occupy singleparticle states is written as: where f (ǫ i ) is the occupation of the eigenenergy ǫ i , where i denotes an eigenstate, the sum is over one-electronic eigenstates, and k B is the Boltzmann constant. In thermal equilibrium state with respect to an electron subsystem, the occupation f (ǫ i ) can be expressed as the Fermi-Dirac distribution, f (ǫ i ) = (1 + e (ǫi−µ)/kB Te ) −1 , where µ is the chemical potential. To simplify the calculations, only the volume dependence at each T e was considered. The face centered cubic (fcc) structure for primitive cells of Cu at different volumes are calculated in a range of T e between 300 and 25000 K. To analyze the volume dependence of S at high T e , we calculated the band structure and the density of states (DOS) at V 0 , which is the equilibrium volume at T e = 300 K, and at 2V 0 . Both calculations were conducted at T e = 25000 K. The relative error between our calculation value V 0 = 12.23Å 3 and an experimental 29 value V exp = 11.81Å 3 is 3.6 %.
These calculations were performed using xTAPP code, 30 in which the electronic entropy S calculation was implemented. The ultra-soft pseudopotential and the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof exchange-correlation functional 31,32 were used. In the ultra-soft pseudopotential, the 3d 10 and 4s 1 states are treated explicitly as valence states. The electronic structures were calculated with a cutoff energy of 1200 eV for the plane-wave basis and the Brillouin-zone k-point sampling of a Monkhorst-Pack mesh with 12 × 12 × 12 k-points for the fcc primitive cell. The number of bands was 13. In the DOS calculation, only a number of the Brillouin-zone k-point sampling was altered to 16 × 16 × 16 k-points.
The calculation results for F , E and −ST e as a function of V at T e = 300, 15000, 20000 and 25000 K are shown in Fig. 1. Fig. 1(a) shows that the curvatures of electronic-free-energy curves become smaller as the val-ues of T e increase. Between 300 and 20000 K, the minimum points that correspond to equilibrium volume at each T e shift to larger values. Finally, between 20000 and 25000 K, the minimum point vanishes. These results are qualitatively consistent with previous studies for tungsten. 33,34 The present results indicate that if atoms can freely change their interatomic distance, such as atoms near a surface, then they cannot be condensed around T e = 25000 K under the assumption of an isothermal process with respect to T e . The validity of this assumption will be discussed later.
To discuss the main contribution for the disappearance of the minimum point in Fig. 1(a), E and −ST e , which are components of F , are plotted in Figs. 1(b) and (c), respectively. Fig. 1(b) shows that the values of E at high T e are larger than these at low T e in the region of V > V 0 . On the other hand, Fig. 1(c) shows that the values of −ST e at high T e are smaller than these at low T e in the region of V > V 0 . Hence, we find that the main contribution for the disappearance of the electronic free energy minimum originates from −ST e . To discuss the reason for this large benefit of S at large volume, the band structures and the DOS at different volumes of V 0 and 2V 0 are plotted in Fig. 2. From Fig. 2(b) and the definition of S (Eq. (2)), the local DOS near the Fermi energy at a large volume of 2V 0 is larger than those at V 0 and this change increases the value of S. This behavior can be easily understood as a decrease of the hopping energy between atoms due to the increased interatomic distance at the large volume. These consideration for the benefit of S were also suggested in a previous study. 34 We conclude that atoms cannot be condensed around 25000 K due to an increase of the electronic entropy S, and consequently non-thermal ablation occurs. We call this physical mechanism the EED model. It should be noted that this explanation does not require the neutrality breakdown, which is denied by the experimental result 21 in the case of metal ablation. Based on these results, we propose the EED model to explain the nonthermal ablation of metals, not only for copper, because the physical explanation given here is expected to be applicable to all metals. Subsequently, a mathematical model that included the EED model was developed for simulation of the ablation depth, and results of Cu films calculated by employing this model are presented. One of purposes of these calculations is to validate the mathematical model and the EED model. The other is to discuss the contribution of S to the ablation depth. T e and T l were calculated by solving the following two-coupled differential equations for the electron (Eq. (3a)) and the lattice (Eq. (3b)) subsystems, 16 Here, C and κ are the heat capacity and the thermal conductivity, respectively. The e and l indices denote the electron and lattice subsystems, respectively. G is the electron-phonon heat transfer constant and S laser is a source term that describes the energy deposition by the laser pulse. The thermal diffusion along a surface can be neglected because the sum of the laser penetration depth δ = 13 nm 35 and the mean free path of electrons δ mfp = 42 nm 36 of Cu is much smaller than the radius of the irradiated laser spot. Accordingly, the threedimensional Eqs. (3a) and (3b) can be reduced to onedimensional equations. In addition, the melting of Cu is neglected in our calculations. This assumption is expected to be suitable in the case of the low-fluence laser irradiation because molten materials were not detected experimentally 7-9 for these laser condition. To accurately calculate time and space evolution of T e and T l , close attention should be paid to determine these parameters in Eqs. (3a) and (3b). T e dependent heat capacity C e (T e ) is obtained by fitting previous calculation results, 37,38 where C e (T e ) is calculated by taking the derivative of the internal energy E(T e ) with respect to T e . According to the Dulong-Petit law, C l = 3.51 J cm −3 is given by V exp . C l can be assumed to be constant above the Debye temperature T D = 343 K 39 so that this value is a good approximation in our simulation, where T l > T D is almost always satisfied. Based on the Drude model, κ e (T e , T l ) = 1 3 v 2 F C e (T e )τ e (T e , T l ) can be derived. Here, v F = 1.57 × 10 6 m/s 36 and τ e (T e , T l ) are the Fermi velocity and the electron relaxation time, respectively. According to the Fermi liquid theory, τ e (T e , T l ) for electrons with energy near the Fermi energy is approximated as: τ −1 e (T e , T l ) = τ −1 ee (T e )+τ −1 el (T l ) = A e T 2 e +B l T l , where A e and B l are typically assumed to be constant. 40,41 B l = 1.98×10 11 s −1 K −1 was determined so as to reproduce experimental value 42 κ e = 3.99 W cm −1 K −1 at low temperature. A e = 2.22 × 10 6 s −1 K −2 was obtained by a recent first-principles calculation. 43 This value is consistent with the experimental result, 41 in which 6.68×10 5 < A e < 2.89 × 10 6 s −1 K −2 was reported. κ l is often neglected for pure metals because this value is much smaller than that of κ e . G = 1.0 × 10 17 W K −1 m −3 was obtained by first-principles calculation, 44 where it was suggested that the T e dependence of G is small at least below 32000 K. In our calculation, this condition is satisfied for a laser peak fluence F < 1.5 J cm −2 . The source term S laser is assumed to be: 15 where z is a depth spacial coordinate, t is the elapsed time, R is the reflectivity, t p is the laser pulse duration time, t 0 is the delay time of the laser, δ b is the ballistic range of electrons, and β = 4 ln 2. R was recently reported to depend on the number of irradiated pulses n and laser fluence F because of the laser-structured surface and the change of the dielectric constant of the irradiated material. 45 Therefore, to compare the simulation results with the experimental data, 26,27 where over approximately a few tens 26 or one hundred 27 laser pulses are irradiated, the change of reflectivity must be considered. However, it is too difficult to consider these effects without experimental data. Therefore, in our calculations, R n (F ) were determined by fitting the experimental reflectivity 45 as R n (F ) = a n ln F +b n , and simulations were conducted for each reflectivity R n (F ) (n = 10, 50, 100). As results of a n and b n calculations, good R n were obtained, of which the root mean square errors were less than 0.02. The laser penetration depth δ(T e ) was determined by a critical point model 46,47 and parameters were obtained from a previous study. 48 The ballistic range was approximated as: 15 δ b (T e , T l ) = τ e (T e , T l )v F . To solve Eqs. (3a) and (3b), the finite-difference methods were used, where time ∆t, and space step ∆z, were 10 as and 1 nm, respectively. The Neumann boundary condition was used. The thickness of the calculated Cu film was 1 µm. The Cu film was irradiated on the front surface by a laser with the laser pulse duration time t p of 100 fs and a wavelength of 800 nm. The delay time t 0 was 4t p and both the initial T e and T l were set to 300 K. It should be noted that the fitting parameters were not used in the simulation.
It was assumed that to cause ablation, the lattice energy E l (T l (t, z)) at a grid of z, which is defined as E l (T l (t, z)) = C l T l (t, z), must overcome the largest values of electronic free energy F coh (T e (t, z)) (Fig. 3) in the region of V > V 0 . At low T e , F coh (T e ) corresponds to the cohesive energy E coh = 47.76 kJ cm −3 , the value of which is the result of our calculation (Fig. 3(b)), and which agrees well with the experimental value 29,49 47.34 kJ cm −3 . In addition to this assumption, we assume that ablation occurs only at a grid of the surface z surf because the bulk cannot expand freely. Taken together, the criterion for ablation can be expressed as the following inequality: Fig. 1(b) shows that the values of the internal energy E at large volume are larger than those of E at V 0 , even at high T e . Therefore, the absorption of the latent heat E late is required for ablation. In the present simulation, E late is assumed as: If the grid of z surf satisfies Eq. (5), then E late (T l ) begins to be absorbed as the latent heat from the electronic subsystem at the grid of z surf . We consider that the delay time t abs = ∆z/v s is required to cause ablation after a grid becomes the grid of z surf because ablation wouldn't occur until passing through a pressure wave, which is created by previous ablation. The velocity of the pressure wave is assumed to be the velocity of  sound, v s = 4760 m s −1 . 50 To represent ablation, the grid of z surf is removed from the simulation, and the grid of z surf +1 becomes the new surface grid after the absorption of E late (T l ).
The calculation results and experimental data 26,27 are plotted in Fig. 4. The thin lines indicate that the dependence on the number of pulses is not large between the 10 and 100th pulses. Moreover, these thin lines indicate that the calculation results with consideration of the electronic entropy S effects are in good agreement with the experimental data 26,27 in the low-laser-fluence region ( < ∼ 5 J cm −2 ), where the effect of non-thermal ablation is expected to be dominant. [7][8][9] However, the gradients of these lines in the high-laser-fluence region (> 5 J cm −2 ) are underestimated. We consider that the disagreement in the high-laser-fluence region is due to a lack of physical mechanics, such as the ejection of liquid droplets by the recoil pressure 7 created by ablation. On the other hand, in the low fluence region, this effect is expected to be little because molten material isn't created [7][8][9] in this region.
The dashed bold line in Fig. 4 represents calculation results when only thermal ablation is considered, in which the effect of S is ignored. In other words, in these calculations, constant cohesive energy F coh = E coh was used instead of F coh (T e ) for the ablation criterion in Eq. (5). This line has no tail in the low-laser-fluence region; therefore, we suggest that non-thermal ablation in the lowlaser-fluence region is caused by the electronic entropy effects. It should be noted that in our simulation, the change of T e at the grid of z surf during t abs is approximately 10-15% due to thermal flux from deeper grids. Therefore, the isothermal process with respect to T e is expected to be conserved at the grid of z surf during t abs .
In summary, the results of FTDFT calculations show the instability of condensed Cu at high T e due to the electronic entropy effect. Furthermore, the results of the band structure and the DOS indicate that the volumedependence in the electron states near the Fermi energy is the main contribution to this instability. Based on these results, we propose the EED model to describe the physical mechanism for the non-thermal ablation of metals. The results of the developed mathematical model calculations, where the effect of the electronic entropy S is included, show that this model can predict the experimental data in the low-laser-fluence region. This strongly suggests that the electronic entropy effect is dominant for non-thermal ablation of metals.

ACKNOWLEDGMENTS
This work was supported in part by the Innovative Center for Coherent Photon Technology (ICCPT) in Japan. Y. T. was supported by the Japan Society for the Promotion of Science through the Program for Leading Graduate Schools (MERIT).