Microscopic mechanism of ultrashort-pulse laser ablation of metals: a molecular dynamics study incorporating electronic entropy effects

TL;DR

This study develops a molecular dynamics model incorporating electronic entropy effects to better understand the microscopic mechanisms of ultrashort laser pulse ablation of metals, revealing high-energy ion emission and sub-nanometer ablation near the threshold.

Contribution

The paper introduces a novel TTM-MD simulation scheme that includes electronic entropy effects, providing new insights into metal ablation mechanisms under ultrashort laser pulses.

Findings

High-energy ion emission observed near ablation threshold

Sub-nanometer depth ablation confirmed by simulations

Electronic entropy significantly influences spallation processes

Abstract

The microscopic mechanism of metal ablation induced by ultrashort laser pulse irradiation is investigated. A two-temperature model scheme combined with molecular dynamics (TTM-MD) is developed to incorporate electronic entropy effects into the simulation of metal ablation while satisfying the energy conservation law. Simulation with the TTM-MD scheme reveals that ultrashort laser pulse irradiation near the ablation threshold causes high-energy ion emission and sub-nanometer depth ablation, as observed experimentally, due to the electronic entropy effect. It is also shown that the electronic entropy effect is also significant in spallation.

Equations

$$C_e^n\frac{d{T}_e^n}{dt}=\nabla \cdot$$

$$-{\sum }_i^{N^n}{\mathbf{v}}_i\frac{\partial }{\partial {\mathbf{r}}_i}\left(S^n{T}_e^n\right)+I^n,$$

$$C_e^n$$

$$G^n$$

$${\kappa }_e^n$$

$$I^n$$

$$N_0$$

$$C_l^n\frac{\text{d}{T}_l^n}{\text{d}t}$$

$$T_l^n$$

$$\frac{d{\mathbf{r}}_i}{dt}$$

$$m\frac{d{\mathbf{v}}_i}{dt}$$

$$F^n=E^n-S^n{T}_e^n.$$

$$\frac{d}{dt}$$

$$=\frac{d{T}_e}{dt}\frac{\partial E}{\partial T_e}+\sum _i^N\frac{d{\mathbf{r}}_i}{dt}\frac{\partial E}{\partial {\mathbf{r}}_i}+\sum _i^N{\mathbf{v}}_i\left(m\frac{d{\mathbf{v}}_i}{dt}\right)$$

$$=\frac{d{T}_e}{dt}\frac{\partial E}{\partial T_e}+\sum _i^N{\mathbf{v}}_i\frac{\partial E}{\partial {\mathbf{r}}_i}-\sum _i^N{\mathbf{v}}_i\left[\frac{\partial (E-S{T}_e)}{\partial {\mathbf{r}}_i}+m\xi {\mathbf{v}}_i\right]$$

$$=C_e\frac{d{T}_e}{dt}+\sum _i^N{\mathbf{v}}_i\frac{\partial }{\partial {\mathbf{r}}_i}\left(S{T}_e\right)-\sum _i^{N}m\xi {\mathbf{v}}_i^2.$$

$$C_e\frac{d{T}_e}{dt}=-\sum _i^N{\mathbf{v}}_i\frac{\partial }{\partial {\mathbf{r}}_i}\left(S{T}_e\right)+\sum _i^{N}m\xi {\mathbf{v}}_i^2.$$

$$\sum _i^{N^n}m{\xi }^n{\mathbf{v}}_i^2+G^n(T_e^n-T_l^n)=0.$$

$$D_{\text{tot}}^n=-\int \nabla \cdot ({\kappa }_e^n\nabla T_e^n)dt,$$

$$I_{\text{tot}}^n=\int I^{n}dt.$$

$$\frac{d}{dt}$$

$$=C_e^n\frac{d{T}_e^n}{dt}+\sum _i^{N^n}{\mathbf{v}}_i\frac{\partial }{\partial {\mathbf{r}}_i}\left(S^n{T}_e^n\right)-\sum _i^{N^n}m{\xi }^n{\mathbf{v}}_i^2$$

$$-\nabla \cdot ({\kappa }_e^n\nabla T_e^n)-I^n$$

$$=C_e^n\frac{d{T}_e^n}{dt}+\sum _i^{N^n}{\mathbf{v}}_i\frac{\partial }{\partial {\mathbf{r}}_i}\left(S^n{T}_e^n\right)+G^n(T_e^n-T_l^n)$$

$$-\nabla \cdot ({\kappa }_e^n\nabla T_e^n)-I^n.$$

$$C_e^n\frac{d{T}_e^n}{dt}=\nabla \cdot ({\kappa }_e^n\nabla T_e^n)-G^n(T_e^n-T_l^n)+I^n.$$

$$E_{\text{cons}}=E_{\text{MD}}+D_{\text{CM}}-I_{\text{tot}}.$$

$$F_{\text{uncons}}=E_{\text{cons}}-\sum _n^{\text{MDcells}}{S}^n{T}_e^n.$$

$$p^n$$

$$\times$$

Full text

I introduction

Laser ablation is widely employed in industry as a method for laser processing (cutting, drilling), pulsed laser deposition, Watanabe, Takahashi, and Tsushima (1998); Yoshitake, Nagamoto, and Nagayama (2000) and nanoparticle production. Fojtik and Henglein (1993); Neddersen, Chumanov, and Cotton (1993) The physical mechanism of laser ablation, especially with ultrashort-pulse lasers (fs laser), has attracted attention in science and industry Kobayashi, Takahashi, and Tani (2020); Kobayashi et al. (2021) because it involves remarkable phenomena that cannot be observed with long-pulse lasers, such as almost no thermal damage, Chichkov et al. (1996); Shaheen, Gagnon, and Fryer (2013) depth of less than 1 nm, Hashida et al. (1999, 2002); Miyasaka et al. (2012) and emission of high-energy ions. Miyasaka et al. (2012); Hashida et al. (2010); Dachraoui and Husinsky (2006); Dachraoui, Husinsky, and Betz (2006) This ablation, which cannot be explained under the assumption of thermal equilibrium, is referred to as non-thermal ablation, whose effects have been reported to be dominant near the ablation threshold fluence. Hashida et al. (2002); Miyasaka et al. (2012); Hashida et al. (2010); Chichkov et al. (1996); Shaheen, Gagnon, and Fryer (2013); Momma et al. (1996) Although tremendous efforts using both experimental and theoretical approaches have been devoted to elucidating the physical mechanism of the non-thermal ablation of metals, discrepancies exist between experiments and previous theoretical simulations.

Molecular dynamics (MD) simulation is a powerful computational tool to elucidate the microscopic mechanisms of metal ablation. To date, MD simulation of laser ablation in the low-laser-fluence region has been reported for several metals [aluminum (Al), Wu and Zhigilei (2013) silver (Ag), Ji and Zhang (2017) copper (Cu), Foumani and Niknam (2018); Schäfer, Urbassek, and Zhigilei (2002) gold (Au), Ivanov and Zhigilei (2003) nickel (Ni), Ivanov and Zhigilei (2003) and platinum (Pt) Rouleau et al. (2014)]. These calculation results proposed the following explanation for the ablation mechanisms of metals caused by irradiation with ultrashort laser pulse with low laser fluence. With irradiation by ultrashort laser pulse near the ablation threshold fluence, the laser-deposited energy raises the surface temperature so that the surface starts to expand and begins to melt. Subsequently, tensile stress occurs near the surface region, and as a result, a molten surface layer is spalled, whose thickness is more than $10\,\text{nm}$. Zhigilei, Lin, and Ivanov (2009); Ji and Zhang (2017); Foumani and Niknam (2018); Schäfer, Urbassek, and Zhigilei (2002); Gan and Chen (2009); Ivanov and Zhigilei (2003) This ablation process is called spallation, and has been observed in experiments. von der Linde, Sokolowski-Tinten, and Bialkowski (1997); Sokolowski-Tinten et al. (1998a, b); Sokolowski-Tinten and von der Linde (2000); Rouleau et al. (2014) As the laser fluence increases, the thickness of the spalled layer decreases, and eventually small clusters and atoms are emitted from the overheated surface. Wu and Zhigilei (2013) This ablation process is called phase explosion, and the main cause of this process is considered to be the thermodynamic instability of the overheated surface. Miotello and Kelly (1999); Bulgakova and Bulgakov (2001)

Previous MD simulation studies have argued that isolated atoms are not emitted with irradiation by ultrashort laser pulse near the ablation threshold, and that the phenomenon that occurs near the ablation threshold is spallation. This means that these explanations have a fatal problem in describing the metal ablation induced by ultrashort laser pulses near the ablation threshold fluence, since the emission of high-energy ions and sub-nanometer depth ablation have been experimentally observed in this fluence region. Hashida et al. (1999, 2002); Miyasaka et al. (2012); Hashida et al. (2010); Dachraoui and Husinsky (2006); Dachraoui, Husinsky, and Betz (2006)

It is considered that this discrepancy between MD simulations and experiments comes from a lack of physical mechanisms in previous MD simulations, where the force acting on atoms is assumed to not be changed even in a highly excited laser-irradiated system. Based on this consideration, some physical mechanisms have been proposed to explain the process of the non-thermal ablation of metals. Tao and Wu (2014); Li et al. (2015); Norman, Starikov, and Stegailov (2012); Norman et al. (2013); Stegailov and Zhilyaev (2015, 2016); Ilnitsky et al. (2016); Miyasaka et al. (2012); Dachraoui and Husinsky (2006); Dachraoui, Husinsky, and Betz (2006) One of the most famous ones is the Coulomb explosion (CE), which has been experimentally verified in the case of a semiconductor Zhao and Shin (2013) an insulator, Stoian et al. (2000a, b) and a molecular system. Sato et al. (2008) CE describes the physical mechanism of non-thermal ablation as follows. Under intense laser irradiation, electrons are emitted from a laser-irradiated surface due to the photoelectric effect and/or the thermionic emission process so that strong Coulomb interaction occurs between positively charged ions at the ionized surface. Hence, when the Coulomb interaction is strong enough to overcome the bonding forces between these ions, they are emitted from the surface. If CE plays a dominant role in the laser ablation process, the peak velocity of the emitted ions is scaled by the valence of the emitted ions, which has been observed by time-of-flight experiments in a semiconductor, Zhao and Shin (2013) insulator, Stoian et al. (2000a) and a molecular system. Sato et al. (2008) These observations have been regarded as conclusive evidence of CE in these materials. On the other hand, the peak velocity of the emitted Cu ions is not scaled by the valence of the ions. Zhao and Shin (2013) In addition, other experimental result Li et al. (2011) showed that the electric field near the surface created by the laser irradiation is shielded within the time duration of the probe pulse ($200\,\text{fs}$), and this fast electrostatic shielding is expected to be natural because the inverse of the plasma frequency is very fast ($<1\,\text{fs}$) in bulk Cu. This confirms the consideration based on the continuum model (CM) calculation Lin et al. (2012) that the electric field near surface due to the electron emission is shielded by high-mobility electrons in the bulk metal before the CE can occur. Hence, the validity of the CE in metals is questionable. Besides CE, other possible origins of the non-thermal ablation of metals have been proposed, for example, the kinetic energy of free electrons and changes in the charge distribution. Norman, Starikov, and Stegailov (2012); Norman et al. (2013); Stegailov and Zhilyaev (2015, 2016); Ilnitsky et al. (2016) However, the validity of these explanations is still under debate.

Recently, we have shown by finite-temperature density functional theory (FTDFT) calculations that the laser-irradiated bulk metal becomes unstable due to the electronic entropy effect. Tanaka and Tsuneyuki (2018) This result suggests that the non-thermal ablation of metals is induced by the electronic entropy effect, and this explanation for the non-thermal ablation of metals is called the electronic entropy-driven (EED) mechanism. Tanaka and Tsuneyuki (2018) Based on the EED mechanism, we have developed a CM where the well-known two-temperature model (TTM) Anisimov, Kapeliovich, and Perel’man (1974) and the electronic entropy effect are incorporated, and succeeded in quantitatively describing the experimental ablation depth Colombier et al. (2005); Byskov-Nielsen et al. (2010) in the low-laser-fluence region.

To further discuss the validity of the EED mechanism and investigate the effect of electronic entropy on the non-thermal ablation of metals, MD simulation is preferred over CM simulation since it can directly describe atom emission and sub-nanometer scale ablation, which are characteristic of non-thermal ablation. Hashida et al. (1999, 2002); Miyasaka et al. (2012); Hashida et al. (2010); Dachraoui and Husinsky (2006); Dachraoui, Husinsky, and Betz (2006) Previously, to elucidate the microscopic mechanism of laser-irradiated metals, a two-temperature model combined with molecular dynamics (TTM-MD) scheme has been employed. Murphy et al. (2015, 2016); Daraszewicz et al. (2013) However, the previous TTM-MD scheme is not appropriate for a system in which the electronic entropy makes a large contribution, and to our knowledge, there have been no TTM-MD schemes that satisfy the law of conservation of energy in such a system. Therefore, to carry out reliable TTM-MD simulation of the non-thermal ablation of metals, where the electronic entropy effect is proposed to be large, Tanaka and Tsuneyuki (2018) it is necessary to develop a new TTM-MD scheme.

The purpose of this study was twofold. The first was to develop a TTM-MD scheme that satisfies the law of conservation of energy even in a system where electronic entropy effects make a dominant contribution. The other was bridging the discrepancy between experiment and previous theoretical simulations regarding the explanation of the non-thermal ablation of metals by elucidating the effect of electronic entropy in these phenomena.

The outline of this paper is as follows. In Sec. II, the developed TTM-MD scheme and computational details are explained. Owing to this development, the TTM-MD simulation can be performed while satisfying the law of conservation of energy even in a system where electronic entropy effects are large. In Sec. III, it is firstly shown that the law of conservation of energy is satisfied with reasonable accuracy in the developed TTM-MD simulation. Subsequently, calculation results for the ultrashort-pulse laser ablation of a Cu film using the developed TTM-MD simulation are exhibited. Here, the microscopic mechanisms of the metal ablation and the effect of the electronic entropy are investigated. To confirm the validity of the TTM-MD simulation and the EED mechanism, the ablation depth in the TTM-MD simulation is compared with previous calculations Tanaka and Tsuneyuki (2018) and experimental data. Colombier et al. (2005); Byskov-Nielsen et al. (2010) Finally, a brief conclusion is provided in Sec. IV.

II Calculation Methods

II.1 Two-temperature model (TTM)

Fig. 1 represents a schematic image of the two-temperature model (TTM), Anisimov, Kapeliovich, and Perel’man (1974) which has been widely used to describe laser-irradiated systems. Daraszewicz et al. (2013); Ernstorfer et al. (2009); Giret, Gellé, and Arnaud (2011); Norman, Starikov, and Stegailov (2012); Norman et al. (2013); Recoules et al. (2006); Inogamov et al. (2012); Wang et al. (2017); Murphy et al. (2015, 2016) Ultrashort-pulse laser irradiation of a metal surface changes the electron subsystem (ES) from the ground state into excited states by the absorption of single or multiple photons. The ES is thermalized to the Fermi-Dirac distribution with the electron temperature $T_{e}$ via the electron-electron (el.-el.) interaction, of which the scattering time $\tau_{ee}$ is approximately $10\mathchar 45\relax 100\,\text{fs}$ in metals. Mueller and Rethfeld (2013); Brown et al. (2016)

In this time scale, the ES and the lattice subsystem (LS) do not reach equilibrium with each other, so $T_{e}$ is higher than the lattice temperature $T_{l}$. Ordinarily, the maximum $T_{e}$ reaches values more than 10 times higher than the final equilibrium temperature ($T_{e}\approx T_{l}$) since the heat capacity of an electron is very much smaller than that of the lattice. $T_{l}$ begins to increase by energy transfer from the ES via electron-phonon (el.-ph.) scattering, for which the relaxation time $\tau_{el}$ is larger than several picoseconds. Schoenlein et al. (1987); Elsayed-Ali et al. (1987, 1991); Hohlfeld et al. (2000) Therefore, under the assumption of the instantaneous and local thermalization in the ES and the LS, the ultrashort-laser-irradiated metals can be described by $T_{e}>T_{l}$ before $\tau_{el}$. This explanation is the main concept of the TTM. Based on the TTM, many previous studies Daraszewicz et al. (2013); Ernstorfer et al. (2009); Giret, Gellé, and Arnaud (2011); Norman, Starikov, and Stegailov (2012); Norman et al. (2013); Recoules et al. (2006); Inogamov et al. (2012); Wang et al. (2017) have been successful in description of the experimental data.

II.2 Two-temperature model combined with molecular dynamics (TTM-MD) scheme

Here, we explain a newly developed calculation scheme for simulating the atom dynamics of metal ablation caused by irradiation with an ultrashort-pulse laser. In the scheme, the MD scheme is hybridized with the TTM scheme to express the non-equilibrium state between the ES and the LS. To decrease the computational cost of large-scale atomistic simulations, the CM is partly employed in LS as well as in ES. Fig. 2 represents a schematic image of this calculation scheme. Hereafter in this paper, this calculation scheme is called the TTM-MD scheme.

In the TTM-MD scheme, electronic effects, such as a highly excited ES near the surface, electronic thermal diffusion, electron-phonon scattering, and energy absorption due to the electronic entropy effect, are incorporated into the MD simulation through the TTM. In the TTM-MD simulation, atom dynamics are calculated based on the MD scheme, and at the same time, other time developments such as that of $T_{e}$ are calculated by employing the TTM. For reduction of computational cost, the CM is also used to calculate the time development of $T_{l}$ deep inside the Cu film (Region 2 in Fig. 2). The dynamics of each atom in this region is not as dominant in the atom dynamics of laser ablation, so only the time development of $T_{e}$ and $T_{l}$ are calculated in this region. This region is called the CM region of LS and plays an important role in the thermal dissipation of the energy deposited by laser irradiation. On the other hand, the region near the surface in which atoms exist is called the MD region of LS (Region 1 in Fig. 2). With the volume change due to expansion or ablation, the position of the surface and the CM region change during simulation. Periodic boundary conditions are used in the $x$-axis and $y$-axis directions in Fig. 2. The free boundary condition is employed between the CM and the MD regions of LS.

The local electron temperature $T_{e}^{n}$ and the local lattice temperature $T_{l}^{n}$ are defined in three-dimensional (3D) cells, where $n$ is the index of the 3D cells. A region surrounded by dotted lines in Fig. 2 represents one of the 3D cells. Although $T_{l}^{n}$ is referred to as the local “lattice” temperature, we do not imply that a crystalline structure is assumed in the 3D cells. In other words, $T_{l}^{n}$ represents not only the lattice temperature but also the temperature of the atoms. Besides, it is noted that $T_{l}^{n}$ of the MD region of LS represents the instantaneous temperature of atoms.

The time development of $T_{e}^{n}$ is calculated by solving the following nonlinear differential equation:

[TABLE]

where $\bm{r}_{i}$ and $\bm{v}_{i}$ are the position and the velocity of atom $i$ in the $n$-th cell, respectively. $C_{e}^{n}$ is the electronic heat capacity, $\kappa_{e}^{n}$ is the electronic thermal conductivity, $G^{n}$ is the electron-phonon heat transfer constant, $N^{n}$ is the number of atoms, $S^{n}$ is the electronic entropy, and $I^{n}$ is the energy deposited by laser irradiation at each $n$-th 3D cell. These quantities are calculated at each 3D cell by the following equations:

[TABLE]

Here, $\rho_{0}$ is the bulk density in the equilibrium states and $V_{c}$ is the volume of each 3D cell. $C_{e}(T_{e}^{n})$, $G$, $\kappa_{e}(T_{e}^{n},T_{l}^{n})$, and $I(T_{e}^{n},T_{l}^{n})$ represent each physical property per unit volume. The values for these properties are the same as those used in the previous study, Tanaka and Tsuneyuki (2018) whose details are explained in the Supplemental Material. The third term on the right-hand side of Eq. (1) represents the absorption of energy by the electronic entropy. The derivation of Eq. (1) based on the law of conservation of energy is explained in Sec. II.3.

$T_{l}^{n}$ in the CM and the MD regions of LS is calculated by solving Eqs. (3a) and (3b), respectively:

[TABLE]

Here, $k_{B}$ is the Boltzmann constant, $\bm{v}_{c}^{n}$ is the average velocity of atoms (center-of-mass velocity) in the $n$-th 3D cell, while $C_{l}^{n}=\frac{N^{n}}{N_{0}}C_{l}$ is the lattice heat capacity in the $n$-th 3D cell. $C_{l}$ is the lattice heat capacity per unit volume, and details are also explained in Supplemental Material.

The atomic dynamics in the MD region of LS are calculated by solving the following equations:

[TABLE]

Here, $m$ is the mass of an atom and ${\xi}^{n}$ is a coefficient that represents the force deriving from the electron-phonon interaction. $F^{n}$ of Eq. (4b) is the free energy of the ES in the $n$-th 3D cell, and the definition is given in the following equation:

[TABLE]

Here, $E^{n}$ represents the internal energy. In this study, $F^{n}$ and $E^{n}$ are calculated using the $T_{e}$-dependent inter-atomic potential (IAP), which is based on the embedded atom method (EAM) potential. The functional form and parameter values for the $T_{e}$-dependent IAP of Cu were proposed in a previous study. Tanaka and Tsuneyuki (2022) The previous study reported that this $T_{e}$-dependent IAP can reproduce the FTDFT results of $T_{e}$-dependent physical properties, such as the volume dependence of $F_{n}$ and $E_{n}$, and the phonon dispersion. Moreover, MD simulations using the $T_{e}$-dependent IAP quantitatively reproduce the results of MD simulation using FTDFT; for example, the time development of the elastic properties of nano-scale slabs and an ablation threshold $T_{e}$.

II.3 Law of conservation of energy

In this study, the developed TTM-MD simulations were performed to investigate the microscopic mechanism of metal ablation induced by irradiation with an ultrashort-pulse laser. Although the $T_{e}$-dependent IAP is used in some previous TTM-MD simulations, Daraszewicz et al. (2013); Murphy et al. (2015, 2016) the law of conservation of energy is considered not to be satisfied, the reason for which is explained below. In these simulations, because the laser has too small a fluence to cause ablation, the effect of electronic entropy may not be very large and deviation from the law of conservation of energy might be negligible. On the other hand, the electronic entropy effect is proposed to be dominant in metal ablation induced by irradiation with an ultrashort-pulse laser. Tanaka and Tsuneyuki (2018) In this case we must carefully take the electron entropy effect into account to realize energy conservation in the TTM-MD simulation. In the following Sec.II.3, it is shown that the law of conservation of energy is satisfied in the developed TTM-MD scheme, theoretically.

First, to simplify the situation, the laser-deposited energy and the energy flow among the 3D cells are neglected. In other words, only energy exchange between the ES and the LS in the 3D cells is considered. In this situation, the conserved energy of the 3D cells is the internal energy: $E^{n}+\sum_{i}^{N^{n}}\frac{1}{2}m\bm{v}_{i}^{2}$.

The time derivative of the conserve energy can be calculated easily as

[TABLE]

Here, to simplify notation, the 3D cell index $n$ is omitted. In the second equality, Eqs. (4a), (4b), and (5) are used. In the third equality, the definition of the electronic heat capacity $C_{e}(T_{e})={\partial E(T_{e})}/{\partial T_{e}}$ is used. Since the time derivative of the conserved quantity is 0, the following equation can be derived:

[TABLE]

The first term on the right-hand side of this equation represents the absorbed energy due to the electronic entropy and the second term is the exchange energy due to the electron-phonon interaction. It is a fundamental assumption in TTM that the electron-phonon interaction is represented by a single linear coupling term of the form of $G^{n}(T_{e}^{n}-T_{l}^{n})$. Hohlfeld et al. (2000) Previously, based on this assumption, the value of $G^{n}$ for Cu has been investigated by experiment Elsayed-Ali et al. (1987) and theoretical calculations, Lin, Zhigilei, and Celli (2008); Migdal et al. (2016, 2016); Petrov, Inogamov, and Migdal (2013); Migdal et al. (2015) for which the details are explained in the Supplemental Material. Therefore, for energy conservation with respect to the electron-phonon interaction between ES and LS, the following equation must be satisfied:

[TABLE]

Subsequently, we added the effect of the electronic thermal diffusion energy $D^{n}_{\text{tot}}$ and laser-deposited energy of the $n$-th 3D cell $I^{n}_{\text{tot}}$ to this scenario. The former effect can be expressed as

[TABLE]

and the latter effect can be written as

[TABLE]

In this situation, the conserved energy in each 3D cell is $E^{n}+\sum_{i}^{N_{n}}\frac{1}{2}mv_{i}^{2}+D^{n}_{\text{tot}}-I^{n}_{\text{tot}}$. Therefore, the time derivative of the conserved energy can be written as

[TABLE]

Here, Eqs. (6), (9), and (10) are used in the first equality, and Eq. (8) is used in the second equality. Since the time derivative of the conserved quantity is 0, Eq. (1) can be derived using Eq. (11). Consequently, Eq. (1) is derived based on the law of conservation of energy.

In previous studies, Daraszewicz et al. (2013); Murphy et al. (2015, 2016) forces acting on the atoms were calculated by the spatial derivative of the free energy calculated using the $T_{e}$-dependent IAP, and the energy exchange due to the electron-phonon interaction was considered. However, the absorbed energy due to the electronic entropy effect was ignored. This means that, in conventional simulations, the time development of $T_{e}$ is calculated by the following equation:

[TABLE]

Hence, energy that is used to accelerate atoms and to raise the internal energy surface is supplied from the virtual electron thermal bath [Fig. 3], because the third term on the right-hand side of Eq. (1) is ignored in the conventional TTM-MD scheme [Eq. (12)]. In this study, we developed the TTM-MD scheme by adding the $-{\sum}_{i}^{N^{n}}\bm{v}_{i}\frac{{\partial}}{{\partial}\bm{r}_{i}}\left(S^{n}T_{e}^{n}\right)$ term to the equation of the conventional TTM-MD scheme, which enabled us to perform simulations that satisfy the law of conservation of energy even in a system where the electronic entropy effects are large.

II.4 Calculation Conditions

Here, we explain the calculation conditions. The lateral dimensions of a laser-irradiated Cu film are $3.615\,\text{nm}\times 3.615\,\text{nm}$, which is ten times the lattice constant of the conventional unit cell of the face-centered cubic (fcc) structure of Cu. The initial MD and CM regions of LS are about $361.5\,\text{nm}$ and $638.5\,$nm, respectively. Hence, the thickness of the computational Cu film is $1\,\mu$m. The total number of atoms in the computational cell is about $4.0\times 10^{5}$. The surface of the film is a (001) free surface of the fcc structure. The laser pulse shape is assumed to be Gaussian. The pulse duration times of an ultrashort-pulse laser and a ps-pulse laser are $100\,\text{fs}$ and $200\,\text{ps}$, respectively. The size of the 3D cells is $1.205\,\text{nm}\times 1.205\,\text{nm}\times 1.205\,\text{nm}$. Therefore, the space step $\Delta x_{\text{CM}}$ is $1.205\,\text{nm}$.

Eqs. (1) and (3a) are solved by a finite difference method (FDM). To solve Eqs. (4a) and (4b), the velocity Verlet algorithm is used. A value of the time step ${\Delta}t$ is $10\,\text{as}$. This value is much shorter than the time step for ordinary MD simulations. To reduce the calculation cost, the time step for MD calculation ${\Delta t}_{\text{MD}}$ is set to ${\Delta t}_{\text{MD}}=n_{\text{MD}}{\Delta}t$, where $n_{\text{MD}}$ is an integer number. In Sec. III.1, we determine a suitable time step $\Delta t_{\text{MD}}$ so that the law of conservation of energy is satisfied with little error.

In our simulations, before irradiation of a laser pulse on the Cu film, the computational cell was relaxed using the Nos$\acute{\text{e}}$-Hoover thermostat Hoover (1985) at $300\,\text{K}$ for $800\,\text{ps}$, where $\Delta t_{\text{MD}}=5\,\text{fs}$ was used. Details of the calculation flow of the TTM-MD are explained in the Supplemental Material.

III Results and Discussion

III.1 Conservation of energy

Here, it is shown that the developed TTM-MD scheme satisfies the law of conservation of energy with a small error. In addition, we show calculation results of the $\Delta t_{\text{MD}}$ dependence of the conserved energy, which is investigated to choose an appropriate $t_{\text{MD}}$ for the following simulations.

Fig. 4 represents a schematic image of the conservation of energy for the MD region ($E_{\text{cons}}$), which is calculated to investigate whether the developed TTM-MD simulation satisfies the law of conservation of energy. $E_{\text{cons}}$ can be written as follows:

[TABLE]

Here, $E_{\text{MD}}$, $D_{\text{CM}}$, and $I_{\text{tot}}$ represent the internal energy of the MD region, the energy thermally diffusing to the CM region, and the energy deposited on the Cu film by the laser, respectively. The internal energy of the MD region is expressed as $\sum_{n}^{\text{MD\,cells}}(E^{n}+\sum_{i}^{N^{n}}\frac{1}{2}m\bm{v}_{i}^{2})$, where the first summation is taken over all MD cells. For comparison, we calculate the free energy of the MD region ($F_{\text{uncons}}$), which is regarded as the conserved energy in the conventional TTM-MD scheme. The definition of $F_{\text{uncons}}$ is as follows:

[TABLE]

Fig. 5 represents the $\Delta t_{\text{MD}}$ dependence of $E_{\text{cons}}$ (solid line) and $F_{\text{uncons}}$ (dotted line) in TTM-MD simulations. Black, green, and blue lines represent calculation results using $t_{\text{MD}}=0.5,1.0$, and $5.0\,\text{fs}$, respectively. Fig. 5(a) represents results for the TTM-MD simulation for $30\,\text{ps}$ without laser irradiation. Figs. 5(b) and (c) represent the results of TTM-MD simulations where the Cu film is irradiated by the ultrashort-pulse laser of (b) $J_{0}=0.4\,\text{J}\,\text{cm}^{-2}$ and by the laser of (c) $J_{0}=0.6\,\text{J}\,\text{cm}^{-2}$. Here, $J_{0}$ represents the laser fluence. Ablation does not occur in (b) the former case; on the other hand, ablation occurs in (c) the latter case. In these two simulations, $T_{e}$ near the surface increases to approximately $20,000\,\text{K}$.

Fig. 5 shows that when sufficiently small $\Delta t_{\text{MD}}$ is used, our simulations satisfy the law of conservation of energy with error of several $10\,\text{meV}\,\text{atom}^{-1}$, and that $F_{\text{uncons}}$ is not conserved. Since Figs. 5(b) and (c) show that $E_{\text{cons}}$ returns back to the initial value at $t>20\,\text{ps}$, where low electron temperatures ($T_{e}\simeq 1,000\,\text{K}$) are realized, these errors of $E_{\text{cons}}$ exist only at high $T_{e}$. The previous study Tanaka and Tsuneyuki (2018) showed that the electronic heat capacity ($C_{e}$) calculated using the $T_{e}$-dependent IAP is slightly overestimated compared to that calculated by FTDFT. In this study, $E^{n}$ is calculated using the $T_{e}$-dependent IAP; on the other hand, the time development of $T_{e}^{n}$ is calculated according to Eq. (6), where $C_{e}$ is calculated by FTDFT. Therefore, $E^{n}$ calculated by the $T_{e}$-dependent IAP is expected to be over-estimated, which can be shown in the calculation results of the laser-irradiated system [Figs. 5(b) and (c)]. Hence, to decrease the error for the law of conservation of energy, it is necessary to develop $T_{e}$-dependent IAP that can reproduce the electronic specific heat of FTDFT with higher accuracy.

Fig. 5(a) shows that $E_{\text{cons}}$ is conserved with little or no error $\Delta t_{\text{MD}}$. On the other hand, Figs. 5(b) and (c) represent that small $\Delta t_{\text{MD}}$ is needed to conserve $E_{\text{cons}}$ when the ultrashort-pulse laser is applied. The reason that energy conservation is not satisfied in the long time step $\Delta t_{\text{MD}}=5.0\,\text{fs}$ can be attributed to high-velocity atoms accelerated by laser irradiation.

In all TTM-MD simulations shown in the following, appropriate $t_{\text{MD}}$ are used after verifying whether $E_{\text{cons}}$ is conserved in each simulation. We carry out the TTM-MD simulations using $\Delta t_{\text{MD}}=1.0\,\text{fs}$ when the irradiation laser fluence is $J_{0}<0.9\,\text{J}\,\text{cm}^{-2}$, and $\Delta t_{\text{MD}}=0.5\,\text{fs}$ is used in the TTM-MD simulations when the irradiation laser fluence is $J_{0}\geq 0.9\,\text{J}\,\text{cm}^{-2}$.

III.2 Microscopic mechanisms of metal ablation

In this section, results and analyses of the TTM-MD simulations of the ultrashort-pulse laser ablation are described. The computational cell for the TTM-MD simulations is shown in Fig. 6. The laser comes from the left side of the Cu film, which consists of MD and CM regions (see Figs. 2 and 6). Fig. 6 and snapshots of the atomic configurations are visualized using Open Visualization Tool Stukowski (2010) (OVITO). The calculation results shown in this section were obtained under conditions in which the fluence of the applied laser changed from $0.54$ to $1.00\,\text{J}\,\text{cm}^{-2}$ while its pulse width was fixed at $100\,\text{fs}$.

III.2.1 Ablation near the ablation threshold: emission of atoms

In the TTM-MD simulations using the $T_{e}$-dependent IAP, Tanaka and Tsuneyuki (2022) emission of an atom is observed when the Cu film is irradiated by a laser pulse with $J_{0}=0.55\,\text{J}\,\text{cm}^{-2}$. Whereas, irradiation of the Cu film by a laser with $J_{0}=0.54\,\text{J}\,\text{cm}^{-2}$ does not cause atom emission. From these results, the ablation threshold fluence is estimated to be $J_{0}=0.55\,\text{J}\,\text{cm}^{-2}$, which is about the same as our previous CM simulation results ($J_{0}=0.47\,\text{J}\,\text{cm}^{-2}$). Tanaka and Tsuneyuki (2018) The kinetic energy of the emitted atom is estimated to be $46.5\,\text{eV}$. This excessively high atom energy is consistent with the experimental value (about $30\,\text{eV}$ Hashida et al. (2010)), which is the most probable energy for the Cu+ emitted on irradiation by a laser with the ablation threshold fluence.

As shown in Fig. 7(a), several atoms are emitted from the surface irradiated by a laser with $J_{0}=0.57\,\text{J}\,\text{cm}^{-2}$. These atomic configurations are simulated using the $T_{e}$-dependent IAP, in which the electronic entropy effect is incorporated. TTM-MD simulations with $T_{e}$-independent IAPs were also performed using $300\,\text{K}$ potential parameters at all temperatures. Snapshots of the atomic configurations obtained by this simulation are shown in Fig. 7(b). As can be seen from this figure, laser irradiation at $J_{0}=0.57\,\text{J}\,\text{cm}^{-2}$, which causes atom emission when simulated with the $T_{e}$-dependent IAP, does not cause atom emission when simulated with the $T_{e}$-independent IAP.

Fig. 8 shows the spatial distribution of $T_{e}$ and $T_{l}$. As can be seen from Fig. 8, at $t=0$, $1$, and $2\,\text{ps}$ there is little difference between the simulation results using the $T_{e}$-dependent IAP and that using the $T_{e}$-independent IAP. Therefore, the atom emission cannot be explained by the thermalized kinetic energy of the atoms. Since previous studies Tanaka and Tsuneyuki (2018, 2022) show that the internal energy ($E$) becomes more attractive at high $T_{e}$, it is considered that the origin of the atom emission comes from the electronic entropy ($-ST_{e}$) effect, which is reported to induce large repulsion forces at high $T_{e}$.

The contribution of the electronic entropy effect can be investigated more directly by calculating the energy absorbed due to the electronic entropy effect. Fig. 9 represents the energy absorbed due to the electronic entropy term of Eq. (5) at each depth. The solid, dashed, and dotted lines in this figure represent results at $t=0.5,1.0,$ and $1.5\,\text{ps}$, respectively. From this figure, it can been seen that the electronic entropy effect is large near the surface.

Furthermore, the electronic entropy effect regarding the pressure is investigated since the laser-induced pressure is considered to be important for the occurrence of spallation and phase explosion. Wu and Zhigilei (2013) Fig. 10 represents the distribution of the local pressure along the $z$ direction ($p_{z}$). According to the simple deviation by Basinski $et\,al$. Basinski, Duesbery, and Taylor (1971) based on the virial theorem, the local pressure $p^{n}$ in the $n$-th 3D cell also can be calculated from the following equation:

[TABLE]

Here, $V^{n}$ and $N^{\text{tot}}$ are the volume of the $n$-th 3D cell and the total number of atoms, respectively. $r_{ij}$ and $f_{ij}$ are the distance and force between the $i$-th and the $j$-th atoms, respectively. The bracket means the time average. In our calculation, the value of $p^{n}$ is averaged within $100\,\text{fs}$. We focus only on the local pressure along the $z$ direction, which is the most important for the ablation dynamics.

Fig. 10(b) shows that a pressure of less than $5\,\text{GPa}$ is created in the simulation using the $T_{e}$-independent IAP. On the other hand, Fig. 10(a) shows that a large pressure ($\sim 35\,\text{GPa}$) is created near the surface in the simulation using the $T_{e}$-dependent IAP, and more than $5\,\text{GPa}$ pressure is created even deep inside the region ($>200\,\text{nm}$). The compressive pressure wave created by atom emission and that created by expansion near the surface would not reach the deep interior region because the velocity of sound of bulk Cu is $4.76\,\text{nm}\,\text{ps}^{-1}$. Linde (2004) Therefore, this high compressive pressure in the deep interior region is considered to arise from the repulsion force between atoms due to the electronic entropy effect.

As shown in Fig. 10(a), large negative pressure (tensile stress), which has the potential to induce spallation, is created near the surface ($\sim 5\,\text{nm}$) in the simulation using the $T_{e}$-dependent IAP. Although a negative pressure is created, laser irradiation with $J_{0}=0.57\,\text{J}\,\text{cm}^{-2}$ does not cause spallation even after the compressive pressure wave reaches the MD/CM boundary ($361.5\,\text{nm}$). Therefore, we thought that a higher negative pressure and a higher temperature were necessary to cause spallation, and we found that ion emission occurred at a lower laser irradiation fluence than the fluence that caused spallation. From the number of emitted atoms, the ablation depth is estimated to be $0.65\,\text{nm}$ for laser-irradiation with $J_{0}=0.57\,\text{J}\,\text{cm}^{-2}$. Not only the ablation at $J_{0}=0.55\,\text{J}\,\text{cm}^{-2}$ but also this result has the potential to explain the non-thermal ablation of metals where sub-nanometer ablation was observed. Hashida et al. (1999, 2002); Miyasaka et al. (2012)

III.2.2 Ablation with a laser fluence slightly larger than the ablation threshold: spallation

Here, we exhibit the results of simulations in which the fluence of the applied laser ($J_{0}=0.7\,\text{J}\,\text{cm}^{-2}$) is a little higher than the ablation threshold ($J_{0}=0.55\,\text{J}\,\text{cm}^{-2}$).

Fig. 11 shows snapshots of the atomic configurations of the Cu film irradiated by the ultrashort-pulse laser. As shown in Fig. 11(b), when the $T_{e}$-independent IAP is used in the simulation, ablation does not occur. On the other hand, Fig. 11(a) shows that ablation occurs when the $T_{e}$-dependent IAP is used. In addition, Fig. 12 shows that the absorption due to electronic entropy is more than $4.0\,\text{eV}\,\text{atom}^{-1}$, which is larger than that of simulation at a lower laser irradiation and cohesive energy. According to the calculation results, the electronic entropy plays an important role in causing atom emission even at this laser fluence.

Fig. 13 shows snapshots of the atomic configurations of [math] to $15\,\text{ps}$ after laser irradiation. As shown in these figures, spallation is also observed for laser irradiation with $J_{0}=0.7\,\text{J}\,\text{cm}^{-2}$. Part of the atomic configuration (depth: $\sim 10\,\text{nm}$) of Fig. 13 is shown in Fig. 14. Since it is thought that the trigger for spallation is tensile stress, the time development of local pressure in the TTM-MD simulation is calculated to investigate the electronic entropy contribution. The space distribution of the local pressure along the $z$ direction is shown in Fig. 15(a). At least within $t=5\,\text{ps}$, the pressure wave passes through a point indicated by the red arrow in Fig. 15(a), at which the surface layer is spalled, and a large negative pressure is created. Owing to the negative pressure, a void begins to be formed around $t=9\,\text{ps}$, and as a result, spallation occurs. In the simulation with the $T_{e}$-independent IAP, spallation does not occur at least within $100\,\text{ps}$, which is enough time for the recoil pressure created near the surface to reach the MD/CM boundary. Fig. 15(b) shows that the negative pressure for the simulation using the $T_{e}$-independent IAP is smaller than that using the $T_{e}$-dependent IAP by one order of magnitude. Since it has been widely accepted that large negative pressure is the origin of spallation, Wu and Zhigilei (2013) we consider that one of the reasons that spallation is not caused in the $T_{e}$-independent IAP simulation is the small negative pressure. The reason for the small negative pressure is considered to be the lack of atom emission and the small internal pressure due to neglecting the effect of electronic entropy. From these results, we conclude that the effect of electronic entropy enhances not only atom emission but also spallation.

III.2.3 Ablation a little higher than the ablation threshold: transition to phase explosion

Here, we exhibit the results of a simulation in which the fluence of the applied laser is $J_{0}=1.0\,\text{J}\,\text{cm}^{-2}$.

Fig. 16 shows snapshots of the atomic configurations of the Cu film irradiated by the pulse laser. In the TTM-MD simulation, the $T_{e}$-dependent IAP is used. The ablation depth of this simulation is $37.8\,\text{nm}$, which is estimated from the number of atoms emitted before the pressure wave reaches the MD/CM boundary. Fig. 16 shows that homogeneous evaporation (see Fig. 16 at $t=10\,\text{ps}$) is observed near the laser-irradiated surface, which is considered to be an indication of phase explosion. These results indicate that as the laser fluence becomes larger, the ablation process transforms from spallation to phase explosion. This result is qualitatively consistent with a previous MD simulation Wu and Zhigilei (2013) in which the $T_{e}$-independent IAP is used.

III.3 Ablation depth

Here, we show the simulation results for the ablation depth and comparison of the results with a previous CM simulation Tanaka and Tsuneyuki (2018) and experimental results. Colombier et al. (2005); Byskov-Nielsen et al. (2010) The ablation depth is estimated from the number of atoms emitted before the pressure wave reaches the MD/CM boundary. Our calculation results for the ablation depth are plotted in Fig. 17. The open circles represent the average ablation depth results from three TTM-MD simulations with different initial thermalization times. The averaged values for the ablation depth at $J_{0}=0.55$ and $0.57\,\text{J}\,\text{cm}^{-2}$ in the TTM-MD simulation are $0.01$ and $0.30\,\text{nm}$, respectively. Since spallation occurs sometimes and does not occur at other times at $J_{0}=0.60\,\text{J}\,\text{cm}^{-2}$, there is a wide range of ablation depth at this laser fluence, as shown by the bars accompanying the circles. Therefore, the ablation depth changes by more than two orders of magnitude around $J_{0}=0.60\,\text{J}\text{cm}^{-2}$. The TTM-MD simulation is qualitatively consistent with experiment, Hashida et al. (1999, 2010) where similar large changes in the ablation depth are observed.

The dotted and dashed lines represent the results of the previous CM calculation Tanaka and Tsuneyuki (2018) including the electronic entropy effect and ignoring the electronic entropy effect, respectively. Triangles and squares represent experimental results. Colombier et al. (2005); Byskov-Nielsen et al. (2010) As shown in Fig. 17, our TTM-MD simulations of ablation depth are in qualitative agreement with previous experimental and calculation studies.

III.4 Pulse-width dependence of ablation threshold

Here, the pulse-width dependence of the ablation threshold is investigated. A previous study Hashida et al. (2002) reported that the ablation threshold fluence of an ultrashort-pulse laser is lower than that of a ps-laser. We investigate whether our simulation can qualitatively reproduce this experimental result.

Figure 18 represents the space distributions of $T_{e}$ and $T_{l}$. The time duration of the ps-laser pulse is $200\,\text{ps}$ and its fluence is $J_{0}=0.70\,\text{J}\,\text{cm}^{-2}$. The peak of the laser pulse reaches the surface at $t=200\,\text{ps}$. In this simulation, ablation does not occur, at least within $300\,\text{ps}$. In the case of ultrashort-pulse laser ($100\,\text{fs}$ laser) irradiation at the same fluence, ablation occurs (see Fig. 11). These calculation results show that the ablation threshold fluence of the ultrashort-pulse laser is lower than that of the ps-laser. This means that the developed TTM-MD simulations can qualitatively reproduce the experimental results of the pulse-width dependence of the ablation threshold.

The reason for the difference between irradiation with the ultrashort-pulse laser and ps-laser can be explained as follows. As can be seen from Fig. 18, in the case of ps-laser irradiation, the difference between $T_{e}$ and $T_{l}$ is small compared with the simulation for ultrashort-pulse laser irradiation (see Fig. 8). In addition, $T_{e}$ reaches only about $3,000\,\text{K}$, which is one order of magnitude lower than with the $100\,\text{fs}$ laser irradiation. At low $T_{e}$, the electronic entropy effect is small, so atom emission and spallation are suppressed compared with ultrashort-pulse laser irradiation. Therefore, ablation does not occur for the irradiation with the ps-pulse laser.

IV Conclusion

The microscopic mechanism of metal ablation induced by irradiation with an ultrashort-pulse laser was investigated.

First, a new TTM-MD scheme was developed considering the electronic entropy effect. To satisfy the law of conservation of energy, the correction term [$-{\sum}_{i}^{N^{n}}\bm{v}_{i}\frac{{\partial}}{{\partial}\bm{r}_{i}}\left(S^{n}T_{e}^{n}\right)$] is added to the conventional equation of the TTM-MD scheme. The energy conservation in the new scheme was verified by simulation of the laser-irradiated Cu film with $T_{e}$-dependent IAP.

With the TTM-MD simulations, laser ablation of Cu films with an ultrashort laser pulse was investigated. The TTM-MD simulation predicts high-energy atom emission and sub-nanometer depth ablation near the ablation threshold fluence ($J_{0}=0.55\,\text{J}\,\text{cm}^{-2}$), which were observed in experiments. This finding bridges the discrepancy between experiments and previous theoretical simulations in explaining the physical mechanism of the non-thermal ablation of metals. Comparing the TTM-MD simulation with the electronic entropy effect and that without this effect, it is found that the electronic entropy plays an important role in atom emission. In the case of the ultrashort laser pulse with $J_{0}=0.7\,\text{J}\,\text{cm}^{-2}$, atom emission and spallation were induced only in the case of the $T_{e}$-dependent IAP, indicating that electronic entropy plays an important role in causing not only atom emission but also spallation. Moreover, the TTM-MD results for ablation depth were in harmony with the CM calculation results and the experimental data, qualitatively. Additionally, the dependence on the pulse width was analyzed. Ablation does not occur with irradiation by 200 ps laser pulse with $J_{0}=0.7\,\text{J}\,\text{cm}^{-2}$ since ps-laser irradiation does not realize excessively high $T_{e}$ in the system. Hence, the ablation threshold fluence of the ultrashort-pulse laser is found to be lower than that of the ps-laser, which is consistent with experiment.

In this paper, using the developed TTM-MD scheme, we demonstrated that the electronic entropy effect plays an important role in the ultrashort-pulse laser ablation of metals, and supports the EED mechanism to explain the non-thermal ablation of metals.

Acknowledgements.

This work was supported in part by the Innovative Center for Coherent Photon Technology (ICCPT) in Japan and by JST COI Grant Number JPMJCE1313 and also by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) Quantum Leap Flagship Program (MEXT Q-LEAP) Grant No. JPMXS0118067246. Y. T. was supported by the Japan Society for the Promotion of Science through the Program for Leading Graduate Schools (MERIT).