{\it Ab initio} approach to the lattice softening of an Al slab driven by collective electronic excitations after ultrashort laser pulse irradiation

TL;DR

This study uses ab initio molecular dynamics to investigate how ultrashort laser pulses affect lattice dynamics in aluminum, revealing non-uniform plasmon screening and changes in phonon behavior.

Contribution

It introduces a novel ab initio approach to analyze laser-induced lattice softening and plasmon effects in aluminum, advancing understanding of non-equilibrium phonon properties.

Findings

Lattice distortion propagation slows with increased laser intensity.

Harmonic force components decrease significantly under laser excitation.

Spatially non-uniform plasmon screening influences lattice dynamics.

Abstract

Recent advances in ultrashort laser pulse techniques have opened up a wide variety of applications in both fundamental physics and industrial fields. In this work, $ab$ $initio$ molecular dynamics simulations based on time-dependent density functional theory revealed a steady deceleration of lattice distortion propagation in an aluminum slab with increasing laser pulse intensity. Analysis of the interatomic force revealed a significant reduction in the harmonic terms and non-monotonic growth of anharmonicity. This behavior was characterized by spatially non-uniform force screening by plasmons, which is missing from Born--Oppenheimer molecular dynamics, and is consistent with the current interpretation of laser-induced periodic structure patterning. This work provides a semi-quantitative criterion for modifying the phonon properties of non-equilibrium systems.

Equations

$$U(F^{(2)},F^{(3)},\mathbf{R})=\sum _{m,n,i,j}{F}_{mn,ij}^{(2)}{\xi }_{mn,i}{\xi }_{mn,j}$$

$$+\sum _{m,n,i,j}{F}_{mn,ijk}^{(3)}{\xi }_{mn,i}{\xi }_{mn,j}{\xi }_{mn,k}.$$

$$V^{\prime }$$

$$\simeq$$

$${\ddot{u}}_{m,z}$$

Full text

Material processing techniques using ultrafast intense laser pulses have been widely used in both fundamental physics and industrial fields Sugioka and Chen (2014). In comparison to processing methods based on nanosecond laser pulses, the use of sub-picosecond laser pulse irradiation reduces the thermal and/or energy diffusion into the surrounding medium, which leads to high energy efficiency and fine spatial resolution during patterning. The realization of greater efficiency and finer resolution than the laser wavelength has been vigorously sought. Although more than three decades have passed since early experimental reports Srinivasan et al. (1987); Küper and Stuke (1987), elucidating the material properties under or after pulse irradiation remains at the cutting edge of condensed matter physics. However, the extremely non-equilibrium and multiscale nature of ablation processes continues to hinder research in this area.

Reducing the interatomic potential and energy diffusion to the medium would be favorable for improving the spatial resolution and energy efficiency during laser patterning. Various mechanisms of lattice property modulation have been proposed to investigate the formation of sub-wavelength structures during ablation processes. It is widely known that such structures are formed and that their spatial periodicity depends on the laser pulse duration. When nanosecond laser pulses are used, the periodicity is close to the incident laser wavelength. This is considered to originate from the interference between the incident and reflected laser light Jain et al. (1981); Keilmann and Bai (1982). Meanwhile, femtosecond laser pulse irradiation generates grating structures whose periodicity is one order of magnitude smaller than the laser wavelength Sakabe et al. (2009); Yasumaru et al. (2003); Borowiec and Haugen (2003); Costache et al. (2003); Reif et al. (2002); Miyaji and Miyazaki (2006); Shimotsuma et al. (2003); Wang et al. (2010). Numerous mechanisms to explain these phenomena have been proposed, some of which have considered the contribution of plasmonic excitations Reif et al. (2002); Borowiec and Haugen (2003); Bhardwaj et al. (2006); Miyaji and Miyazaki (2008); Vorobyev et al. (2007); Bonse et al. (2009). The direct measurement of lattice properties during ultrafast processes is difficult and fundamental quantities such as the interatomic force constants are only given for thermalized equilibrium system Kresh et al. (2008). Consequently, further information regarding laser–matter interactions would be of great value.

Since ablation is a multiscale process, previous theoretical approaches have ranged from macro- to microscopic models and the descriptions of the systems have also varied. Hydrodynamic models with a nanosecond time scale and micrometer spatial scale have been reviewed previously Schultz et al. (2013). Classical molecular dynamics (MD) simulations have also been performed for sub-micrometer-scale structures of metals Ivanov et al. (2013); Norman et al. (2012). Many studies using quantum mechanical approaches have assumed thermalization of the subsystems Brown et al. (2016); Waldecker et al. (2016). In a recent advance, a first-principles study highlighted the importance of electronic enthalpy in ablation processes based on finite-temperature density functional theory (DFT) Tanaka and Tsuneyuki (2018).

In this Letter, we demonstrate the application of ab initio Ehrenfest molecular dynamics (EMD) simulations based on time-dependent density functional theory (TDDFT) to investigate the laser-driven suppression of interatomic forces and the volume expansion of crystals. The investigated system was a thin slab comprising nine layers of Al atoms. The (111) surfaces of the Al fcc structure is exposed to the vacuum. We applied an ultrashort laser pulse of infrared light with a wavelength of 800 nm, a full width at half-maximum (FWHM) of 10 fs, and a field amplitude ranging from 0.0 to 3.0 $\mathrm{V/\AA}$. To analyze the force field, we employed a quasi-one-dimensional model in which atoms were coupled to their first nearest neighbors via a potential expanded by the third order of interatomic distance. The force constants were fitted such that the model reproduced the EMD trajectory. We confirmed the volume expansion of the slab and significant suppression of the harmonic terms, which amounted to 38% of the initial values when averaged over the layers after irradiation. The corresponding anharmonic terms were also determined. The force suppression was significant on the surface layers, in contrast to the case of Born–Oppenheimer MD (BOMD) simulations with finite electron temperature. Based on a phenomenological analysis, we attributed this discrepancy to force screening by plasmonic excitations, which is absent from the BOMD framework.

We first introduce our theoretical approach to treating the ultrafast dynamics of the electronic and lattice systems. The theoretical description of the real-time evolution of the electronic system was based on TDDFT. In TDDFT, the electronic state at each time step is obtained by solving the time-dependent Kohn–Sham equation for one-particle orbitals Runge and Gross (1984). We used the local-density approximation (LDA) with a Perdew–Zunger-type exchange-correlation functional Perdew and Zunger (1981). The fourth-order Suzuki–Trotter-type time evolution operator Suzuki (1992); Sugino and Miyamoto (1999) was used to ensure numerical accuracy and unitarity of time evolution. Potentials between time steps were interpolated using the railway curve interpolation scheme for numerical accuracy and time reversibility Sugino and Miyamoto (1999).

In Fig. 1 we show the crystal structure of the thin aluminum slab used in the present simulation. We took the $xy$ plane as parallel to the slab. Fig. 1(a) depicts the in-plane hexagonal unit-cell structure and Fig. 1(b) shows a cross section of the slab along the $z$ axis. The (111) surfaces of the fcc Al crystal were exposed to vacuum layers and the slab was composed of nine atomic layers. The lengths of the $a$ ($b$) and $c$ axes were 5.303 bohr and 60.55 bohr, respectively. The cell parameters were fixed in the EMD simulations. Nine Al atoms were contained in the unit cell. We used the Troullier–Martins-type pseudopotential Troullier and Martins (1991) and a 16$\times$16$\times$1 $k$-point mesh. The plane-wave vector cutoff was 35 Ry for the basis set and 562.5 Ry for the charge density distribution. The time step for the time-dependent simulation was 3.63 attoseconds. The electronic system was coupled to the external field by the length gauge $V=e\bm{r}\cdot\bm{E}$, where $e$ is the electronic charge, $\bm{r}$ is the electron position, and $\bm{E}$ is the external electric field with the polarization vector parallel to the $z$ axis. We applied an ultrashort laser pulse that can be analytically expressed as the product of a Gaussian and a sinusoidal function. The FWHM was 10 fs and the frequency was 375 THZ, which corresponds to a wavelength of 800 nm. The maximum field amplitude ranged from 0.0 to 3.0 V/$\mathrm{\AA}$ in 0.5 V/$\mathrm{\AA}$ intervals.

We conducted simulations to measure the propagation time for in-plane atomic displacement, which was initially induced to an outermost atom. This quantity was defined as the time required for the in-plane displacement of the opposite outermost atom to show its first peak. The left panel of Fig. 1(d) shows the propagation time for various field amplitudes in the EMD simulations. For comparison, we performed BOMD simulations with Fermi–Dirac smearing using the Quantum ESPRESSO package Giannozzi and et al. (2009). We used the Troullier–Martins-type pseudopotential Troullier and Martins (1991), a 16$\times$16$\times$1$k$-point mesh, a plane-wave cutoff of 20 Ry, a charge density cutoff of 80 Ry, and a time step of 0.48 fs. Provided that the smearing adequately approximate the effect of finite electronic temperature, we conducted our BOMD simulations with different smearing widths from 0.01 to 0.1 Ry, as shown in the right panel of Fig. 1(d). The discrepancy between the EMD and BOMD simulations at the lowest field amplitude and smearing width may be attributable to the different methods used to assign the occupations and parameter settings for the pseudopotentials, although it was outside the scope of this study to resolve this discrepancy by fine-tuning these parameters. The steady increase in the propagation time indicates the reduction of the interatomic potential. For quantitative analysis, we constructed a quasi-one-dimensional model in which each atomic layer was coupled to its first nearest neighbor via an interatomic potential $U$ expanded by the third order of interatomic distance:

[TABLE]

Here, the interatomic distance $\xi_{mn,i}=u_{m,i}-R_{m,i}-u_{n,i}+R_{n,i}$, where $u_{m,i}$ is the $i$-th component of the $m$-th atomic coordinate and $R_{m,i}$ is the corresponding equilibrium position. Fig. 1(c) shows a schematic expression of this model. $F^{(2)}=\{F^{2}_{mn,ij}\}$, $F^{(3)}=\{F^{(3)}_{mn,ijk}\}$, and $\bm{R}=\{R_{m,i}\}$ are fitting parameters. We projected the potential $U$ such that it satisfies hexagonal symmetry. Thus, Eq. (1) is reduced to $U(F^{(2)},F^{(3)},\bm{R})=\sum_{mn}\{F^{(2)}_{mn,xx}(\xi_{mn,x}^{2}+\xi_{mn,y}^{2})+F^{(2)}_{zz}\xi_{mn,z}^{2}+F^{(3)}_{xxy}(\xi_{mn,x}^{2}\xi_{mn,y}-\xi_{mn,y}^{3}/3)+F^{(3)}_{mn,xxz}(\xi_{mn,x}^{2}+\xi_{mn,y}^{2})\xi_{mn,z}+F^{(3)}_{zzz}\xi_{mn,z}^{3}\}$, which is characterized by five independent force constants. The potential $U$ is assumed to be invariant under inversion of the $z$ axis $\xi_{mn,z}\rightarrow-\xi_{mn,z}$. In our fitting procedure, the evaluation function was defined as the square of the difference between the acceleration of atoms extracted from the EMD trajectory and those constructed by $F^{(2)}$, $F^{(3)}$, $\bm{R}$. For the fitting, we ran the EMD simulation for 392 fs and randomly selected 50 snapshots of acceleration to construct the evaluation function. All of the atomic positions were initially displaced from their equilibrium positions by 5% of the lattice constant. The displacement vector was set to be antiparallel to those of neighboring atoms. We omitted the high-frequency component of acceleration from the EMD trajectory by applying a 100 THZ cutoff prior to parameter optimization using the Fletcher–Reeves optimization method. Additional details are provided in Sec. S.I of the Supplemental Material Katow and Miyamoto (2019).

In Fig. 2(a), (b), and (c), we show the fitted second-order force constants and interlayer distances obtained from the EMD simulations. We compared the results with those from the BOMD simulations, as shown in Fig. 2 (d), (e), and (f). The total time for the BOMD simulation was 1.45 ps. We observed significant suppression of $F^{(2)}$ in both cases. The reductions in $F_{zz}^{(2)}$ averaged over the layers amounted to 38% and 56% of the initial values for the EMD and BOMD simulations, respectively. Although there is no direct experimental report for observing these quantities, a neutron scattering experiment involving bulk aluminum revealed a 4.8% reduction in phonon frequency and 10% reduction in force constants for first nearest neighbors when the temperature was increased from 10 to 775 K Kresh et al. (2008). It is counterintuitive that the decrease in $F^{(2)}_{xx}$ appears to saturate while the propagation time of atomic displacement exhibits a steady increase in Fig. 1(d). The absolute values of $F^{(3)}_{xxy}$ and $F^{(3)}_{xxz}$ averaged over the layers are plotted in Fig. 3, revealing steady growth of these quantities with increasing field amplitude after the reduction in $F^{(2)}_{xx}$ became moderate. We can deduce that the delay in propagation was partially due to the increase in $F^{(3)}$ at strong field intensity. The values of $F^{(3)}_{zzz}$ did not converge under our optimization conditions and are therefore not shown. The spatial dependency of $F^{(3)}$ is summarized in Fig. S.1 of Sec. S.II of the Supplemental Material Katow and Miyamoto (2019). Although the spatial non-uniformity of the force constants was large, i.e., the finite size effect was significant, our results provide a semi-quantitative criterion for constructing models in larger systems under extremely non-equilibrium conditions. Strongly enhanced suppression of harmonic terms on surface layers in EMD indicates the emergence of excited electronic states missing from the BOMD framework and we discuss this point next. Hereinafter, we restrict our discussion to the non-uniform force reduction of $F^{(2)}_{zz}$.

We show the frequency spectrum of the Hellmann–Feynman force $|f_{\omega,m}|$ for the $m$-th atom along the $z$ direction in the right panel of Fig. 4(a) for a maximum field amplitude of 3.0 V/$\mathrm{\AA}$. This spectrum was obtained by averaging 1000 spectra of 60 fs long MD data randomly sampled from the last 360 fs of the 392 fs long MD simulation. In the sub-petahertz region, peaks commensurate with the frequency of the laser pulse $\omega_{ph}$ = 375 THZ and its integer multiples $m\omega_{ph}$, namely, the high harmonic oscillation (HHO), up to $m=3$ were confirmed. Furthermore, in the region above 2 PHZ, a very large peak was observed. It is plausible to regard this as plasmonic oscillations relative to the frequency of the volume plasmon and surface plasmon.

We next examined the screening effect of force constants induced by these high-frequency components. It can be easily verified that, when a pair of harmonic oscillators is linearly coupled, the frequency of one is screened while that of the other remains almost unchanged, if the frequency ratio of the two oscillators is very large. We generalize this concept and consider a phenomenological interaction $V^{\prime}$ where the electronic charge density $\rho(\bm{r},t)$ linearly couples to the atomic position $u_{m,z}(t)$ with a coupling constant $\beta_{m}(\omega)$ as

[TABLE]

where $\rho(\bm{r},t)=\int d\omega\rho_{\omega}(\bm{r},t)=\int d\omega\rho_{\omega}(\bm{r})e^{-i\omega t}$ and the index $m$ indicates atoms. We also used the notation $\beta_{m}(\omega)=\beta(\bm{u}_{m},\omega)$. In Eq. (3), we approximated the space integrals of $\beta(\bm{r},\omega)\rho_{\omega}(\bm{r},t)$ by their values at the $m$-th atomic position as $\beta(\bm{u}_{m},\omega)\rho_{\omega}(\bm{u}_{m},t)\Omega$, where $\Omega$ is a fitting parameter. This procedure corresponds to approximating the ion-charge density interaction by a box potential and neglecting the spatial dependency of $\rho_{\omega}(\bm{u}_{m},t)$; thus, $\Omega$ represents the volume of the box potential. Since the frequencies of the HHO and plasmonic peaks in Fig. 4(b) are far higher than typical phonon frequencies, the screening effect can be approximated by $\gamma_{m}(\omega)$ in the equation of motion for the $m$-th atom as follows:

[TABLE]

where $\gamma_{m}(\omega)=\Omega^{2}\beta_{m}(\omega)^{2}/\omega^{2}$, and $\omega_{1}$ and $\omega_{2}$ are frequency cutoffs. We abbreviate the third-order terms as $O(F^{(3)})$. To derive Eq. (4), we assumed that the harmonic-potential-type restoring force $-\omega^{2}\rho_{\omega}$ acts on $\rho_{\omega}$ in its classical equation of motion as $\ddot{\rho}_{\omega}(\bm{r},t)=-\omega^{2}\rho_{\omega}(\bm{r},t)-\partial V^{\prime}/\partial\rho_{\omega}(\bm{r},t)$. The detailed derivation of $\gamma_{m}(\omega)$ is provided in Sec. S.III of the Supplemental Material Katow and Miyamoto (2019). We used the value of the charge density integrated over the $xy$ plane of the unit cell to obtain the linear density per bohr at the $m$-th atomic position. The coupling $\beta_{m}(\omega)$ was computed using $\Omega\beta_{m}(\omega)=<\rho_{\omega}(\bm{u}_{m})^{*}f_{\omega,m}>/<|\rho_{\omega}(\bm{u}_{m})|^{2}>$, where $<\cdots>$ denotes taking the average over randomly sampled spectra as discussed earlier to determine $|f_{\omega,m}|$ in the right panel of Fig. 4(a) obtained by EMD. We show $\Omega\beta_{m}(\omega)$ in the left panel of Fig. 4(a) and $\gamma_{m}=\int d\omega\gamma_{m}(\omega)$ in Fig. 4(b) for $E=$3.0 V/Å by assuming a common $\Omega$ for all atoms. The values for $\gamma_{m}$ at lower field amplitude are shown in Fig. S.2 of Sec. S.III of the Supplemental Material Katow and Miyamoto (2019). We computed $\gamma_{m}$ for $\omega_{1}$ = 300 THZ and 1300 THZ, where the latter case omits the contribution of optical frequency $m\omega_{ph}$. The screening $\gamma_{m}$ was enhanced on the surface atoms, mostly due to the plasmonic component whose peak positions were lower than those of the inner layers as shown in Fig. 4(b). This behavior coincides with the spatial non-uniformity of $F^{(2)}_{zz}$ in Fig. 2. Thus, we conclude that surface-enhanced plasmonic screening of the interatomic force caused the non-uniform spatial dependency. Although we expect that the plasmonic excitation and HHO also contribute to the behavior of $F^{(3)}$, clarifying these effects will require consideration of the higher-order coupling of $u_{m,z}$ and $\rho_{\omega}(\bm{r},t)$ in our model. This will be investigated in our future work.

Thus far our analysis has clarified the significance of plasmonic effects for modeling the evolution of ablation processes at the sub-picosecond time scale. By deducing the physical origin of the harmonic force constant reduction to the plasmonic excitations, we can discuss possible finite size effects in larger systems that are too computationally expensive to treat. Increasing the slab thickness will cause red shift of the plasmonic peaks as it weakens the confinement effect. This may enhance the screening effect of interatomic force owing to the $\omega^{-2}$ dependency of $\gamma_{m}(\omega)$. Weak confinement will also make the spatial dependency of the $F^{(2)}$ reduction rather moderate.

Interatomic force constants are one of the most fundamental quantities of lattice systems, upon which the micro- and macroscopic quantities of crystals, such as the dispersion and lifetime of phonons, heat capacity, and diffusion coefficient of energy, rely. Laser-induced modulation of these quantities is critical to understanding laser ablation processes, and the current work has quantified the modulation of the force constants for both harmonic and anharmonic terms for the first time based on the TDDFT approach. At this ultrafast timescale and non-equilibrium conditions, collective electronic excitations such as plasmons and HHO take the place of the thermalized electrons that play the main role in ordinary BOMD. According to our analysis, the non-uniformity of the interatomic force reduction can be ascribed to the non-uniform force screening by plasmons. This interpretation is consistent with the plasmon-driven mechanism of periodic structure formation at the sub-wavelength scale during ablation processes Reif et al. (2002); Borowiec and Haugen (2003); Bhardwaj et al. (2006); Miyaji and Miyazaki (2008); Vorobyev et al. (2007); Bonse et al. (2009). The investigation of larger systems would be of great interest to us. However, at present such studies are hindered by high computational cost, and hence a phenomenological model may need to be developed to describe the force screening effect. TDDFT is one of the most promising approaches for constructing such models.

This paper is based on the results obtained from the NEDO project “Development of advanced laser processing with intelligence based on high-brightness and high-efficiency laser technologies” (TACMI project). The numerical results described in this Letter were obtained using the supercomputing resources at the Cyberscience Center of Tohoku University.