# Theory of Out-of-Equilibrium Ultrafast Relaxation Dynamics in Metals

###### Abstract

Ultrafast laser excitation of a metal causes correlated, highly nonequilibrium dynamics of electronic and ionic degrees of freedom, which are however only poorly captured by the widely-used two-temperature model. Here we develop an out-of-equilibrium theory that captures the full dynamic evolution of the electronic and phononic populations and provides a microscopic description of the transfer of energy delivered optically into electrons to the lattice. All essential nonequilibrium energy processes, such as electron-phonon and phonon-phonon interactions are taken into account. Moreover, as all required quantities are obtained from first-principles calculations, the model gives an exact description of the relaxation dynamics without the need for fitted parameters. We apply the model to FePt and show that the detailed relaxation is out-of-equilibrium for picoseconds.

Excitation of a material with an intensive, ultrashort optical pulse brings the material’s electrons into a strongly out-of-equilibrium state which is immediately followed by intense, correlated dynamics between the electrons and other degrees of freedom in the material, such as the lattice or spin systems Fujimoto et al. (1984); Beaurepaire et al. (1996); Koopmans et al. (2010); Liao et al. (2016). The ability to measure the ultrafast relaxation dynamics of the involved subsystems using pump-probe techniques has led to the discovery of many unexpected phenomena, such as ultrafast demagnetization Beaurepaire et al. (1996), change of magnetic anisotropy Hansteen et al. (2005), or coherent generation of magnetic precession Ju et al. (1999); van Kampen et al. (2002). More recently, ultrafast generation of lattice strain waves Kim et al. (2012, 2015); Henighan et al. (2016), coherent control of atomic and molecular dynamics Rabitz et al. (2008), coherent phonon generation Harmand et al. (2013); Lindenberg et al. (2017), and laser-induced superconductivity at high temperature Mitrano et al. (2016) have been reported. In spite of the importance of the material’s nonequilibrium state in these laser-induced phenomena, it is surprising that most theoretical descriptions of the ensuing out-of-equilibrium dynamics are based on the widely-used two-temperature model (2TM) Kaganov et al. (1957); Anisimov et al. (1974), which assumes that the electronic and phononic subsystems are separately in thermodynamic equilibrium Allen (1987).

Research on nonequilibrium states of matter has emerged recently as an important area in condensed matter physics (see, e.g., Stojchevska et al. (2014); Aoki et al. (2014)), and consequently several improved models have been developed which incorporate aspects of out-of-equilibrium electronic dynamics Carpene (2006); Mueller and Rethfeld (2013); Carva et al. (2013); Baranov and Kabanov (2014). However, these still lack a complete out-of-equilibrium description of the full system and its time evolution, and contain parameters that are either fitted experimentally or chosen from a macroscopic system at equilibrium. Additionally, recent investigations Henighan et al. (2016); Waldecker et al. (2016, 2017) have emphasized that the assumption of thermal phonons could lead to marked disagreement with experimental observations.

Here we propose a general theory to describe the ultrafast dynamics triggered by ultrashort laser pulses in metals. Contrarily to the two-temperature model, our model employs an out-of-equilibrium description of the electronic and phononic populations and provides the full dynamic description of the nonequilibrium relaxation processes. The decisive quantities that govern the dynamical relaxation are the phonon mode dependent electron-phonon coupling constants and the mode dependent phonon-phonon scattering terms. It is important to stress that both quantities have an explicit dependence on the electronic temperature, the phonon momentum and branch, and, moreover, that they are fully derived from microscopic ab initio theory. Thus, as the theory only uses quantities obtained from ab initio calculations, it provides a parameter free description of the relaxation dynamics and could hence become of great value for modeling and predicting nonequilibrium dynamics following ultrafast laser excitation. As an example, the model is used to describe the ultrafast dynamics in FePt after femtosecond laser irradiation, illustrating as well the limitations of the 2TM.

To describe the nonequilibrium time evolution of the electronic and phononic degrees of freedom of the laser-excited material, we divide the out-of-equilibrium metallic system in different, independent subsystems that interact with one another, schematically shown in Fig. 1. Specifically, we describe the lattice by dividing it into independent phonon subsystems, each of them corresponding to a specific branch and momentum . They interact with one another through phonon-phonon scattering and with the electrons via electron-phonon scattering. These interactions are phonon momentum and branch dependent. Therefore these phonon subsystem populations evolve separately during the nonequilibrium dynamics and we can define a “lattice temperature” (with ) for each of them (the meaning of this definition is discussed further below). The impulsive laser excitation brings a part of the electrons into a nonthermal state (Fig. 1). For the electrons we follow the description of Carpene Carpene (2006), in which the electronic system is divided in a thermal bath that contains the majority of thermal electrons with temperature and in a laser-excited nonthermal electron distribution, which relaxes driving energy into the thermal distribution through the electron-electron and electron-phonon interactions. The latter interaction conducts the energy from the laser-excited electrons to the different lattice subsystems, bringing them to different temperatures, i.e., into a nonequilibrium state. The phonon-phonon scattering causes the transferred energy to be shared between the phonon subsystems, guiding them toward a common lattice temperature, and therefore to a thermal equilibrium of the lattice. Hence, the rates of electron-phonon and phonon-phonon scattering are the quantities that determine the system’s temporal evolution to equilibrium.

To achieve a theoretical formulation we make use of the conservation of total energy and classical kinetic theory. The total lattice energy is given by and the total electronic energy by , where is the phonon energy, the electron Bloch energy not () (with ), and and are phonon and electron populations, respectively. The latter are in local equilibrium given by the Fermi-Dirac and Bose-Einstein distributions

(1) |

We can define the rates of energy exchange as

(2) | |||||

(3) |

where the dot stands for the time derivative and the subscripts denote the different scattering processes that change the distribution. The equivalence in Eq. (3) stems from the conservation of total energy. The time derivatives of the distribution functions due to different scattering terms can be derived from the classical kinetic theory by using the well-known Fermi’s golden rule of scattering theory. By doing so, we obtain an extended version of the Bloch-Boltzmann-Peierls equations (see Ziman (1960), and Supplemental Material (SM) SM ()),

(4) | |||||

(5) | |||||

(6) | |||||

Here and are the electron-phonon and phonon-phonon matrix elements, respectively.

Equations (2) and (3) describe the energy flow between electron and phonon subsystems under the assumption that the diffusion term can be neglected, which is a valid assumption on the short time scale of the typical out-of-equilibrium process. The laser driving field induces the nonequilibrium electronic distribution and enters in the rate equations as source term. Substituting Eqs. (4)–(6) in (2) and (3), and making an expansion of the distribution functions to second order in the phonon mode and electron temperature differences (see SM SM ()), we obtain a set of coupled differential equations which connect the time evolutions of the electron temperature and the temperatures of the different phonon modes,

(7) | |||||

(8) |

and are the temperature-dependent electronic and phonon mode dependent specific heats, respectively, and and are the mode-dependent electron-phonon coupling constants and the mode-dependent phonon linewidths which are caused by the phonon-phonon interactions. and describe the rate transferred from the laser-induced nonequilibrium electron distribution into the thermal electronic bath via electron-electron interaction and into the lattice through electron-phonon interaction, respectively Carpene (2006). is a function of the mode-dependent phonon frequencies and temperatures, which accounts for the second-order term in temperature differences; its full form is given in the SM SM ().

To obtain a full solution of the out-of-equilibrium model, defined by Eqs. (7) and (8), it is necessary to compute the material specific quantities, the phonon and electronic specific heats and phonon mode-dependent electron-phonon and phonon-phonon linewidths. These can be conveniently calculated using the spin-polarized density functional theory (DFT) in the local-density approximation (LDA). Here we have employed the electronic structure code ABINIT Gonze et al. (2009). Specifically, the mode-dependent electron-phonon linewidths were computed as response function within the density functional perturbation theory whereas the mode-dependent phonon-phonon linewidths due to phonon-phonon interaction were determined using many-body perturbation theory in a third-order anharmonic Hamiltonian which includes up to three-phonon scattering Togo et al. (2015). The coupled equations (7) and (8) are then solved numerically, with the ab initio quantities, and without any free fitting parameters.

We emphasize that the phonon branch and wavevector dependent phonon temperatures act here as an auxiliary quantity that allows us to use for each phononic subspace a Bose-Einstein distribution with a local temperature . While this temperature might not be measurable, recent electron diffraction experiments showed that nonequilibrium phonon populations in reciprocal space can be measured Chase et al. (2016). The electron temperature is conversely a quantity that can be obtained from pump-probe photoemission measurements Sun et al. (1993); Guo et al. (2001); Lisowski et al. (2004). As will become evident below, the out-of-equilibrium model [Eqs. (7) and (8)] leads to results that are markedly different from those of the 2TM. Nonetheless, it can be recognized that under simplifying assumptions the 2TM can be obtained from Eqs. (7) and (8). To this end, it is important to realize that is the nonequilibrium equivalent of the electron-phonon coupling constant that is commonly used in the 2TM. Assuming a single phonon temperature for all phonon modes and wavevectors and neglecting quadratic terms in the temperature difference , along with setting , and using (valid under the assumptions of the 2TM, see SM), one recovers the common 2TM Anisimov et al. (1974).

To recognize the importance of the mode and wavevector dependent electron-phonon coupling constant we perform ab initio calculations of this quantity for ferromagnetic FePt, which is the prime material for high-density optic and magnetic recording Lambert et al. (2014); John et al. (2017). Bulk FePt orders in the L1 structure in which the (001) planes are alternatively occupied by Fe and Pt atoms. Our ab initio calculation of the ground state magnetic properties of FePt is in good agreement with previous work FeP ().

In Fig. 2 we show our ab initio calculated mode-dependent electron-phonon coupling constants of FePt at 300 K together with the phonon dispersions along high-symmetry lines in the simple tetragonal Brillouin zone (BZ). We can readily see that the electron-phonon coupling constants are larger for optical phonon modes, reaching several orders of magnitude differences between some points of the BZ (see e.g. optical phonons at the X point compared to acoustic phonons at the point). These findings demonstrate the limitations of considering a single with the lattice at one local thermal equilibrium, since the range of values of is several orders of magnitude. On account of the different coupling strengths the laser-excited electrons will couple mainly to the optical phonon modes.

To account for the scattering processes of electrons away from the Fermi surface, we have additionally included an electron temperature dependence of the mode-dependent electron-phonon coupling constant, using Lin et al. (2008)

(9) |

with being the electron density of states at energy . In Fig. 3 we illustrate the relevance of the electron temperature for . The calculation shows a fast growth of with (red line), as compared with the temperature independent value (yellow line) and an estimated value of used recently (Mendil et al., 2014) to reproduce experimental data with a 2TM (purple line). At high electronic temperatures, is an order of magnitude larger, which will clearly influence the early dynamics of the system after laser irradiation.

Also, we would like to stress that to correctly assess the dependence of the electron heat capacity on the electronic temperature, in our model we calculate it as the derivative of the total electron energy density against the electron temperature Lin et al. (2008), . This is in contrast to the commonly used calculation of from the Sommerfeld expansion of the free energy, which provides a linear temperature dependence, i.e. . The difference between these two different descriptions is shown in the inset of Fig. 3.

Another key quantity of our model is the explicit incorporation and calculation of the phonon anharmonicities, that enter into the model through the mode-dependent phonon linewidths . This quantity is related to the mode-dependent phonon lifetime due to phonon-phonon interaction, via . Notably, in the 2TM it is assumed that such interactions are very strong and lead to very short phonon lifetimes, and therefore to an immediate equilibration of the lattice. In Fig. 4 we show the calculated of FePt at 300 K along high-symmetry lines in the BZ. It can be clearly seen how different phonon modes and branches have different lifetimes, as indicated in the color changing when moving along the phonon dispersions. The phonon lifetimes for the optical branches range from 2 to 10 ps, while the acoustic branches have phonon lifetimes larger than 10 ps, diverging at the point. These timescales are much larger than the initial ultrafast dynamics of the system. Thus, our calculations show that the assumption of an immediately thermalized lattice, as made in the 2TM, is not tenable, since the phonon-phonon lifetimes are at least one order of magnitude larger than the electron-phonon lifetimes.

The combination of the main results shown in Figs. 2 and 4, namely strong mode-dependent excitation of phonon modes via electron-phonon coupling and slow lattice thermalization through phonon-phonon interaction, suggests that the lattice remains out-of-equilibrium not only at sub-ps timescales, but even on much larger timescales. This is a striking difference with respect to the model by Waldecker et al. Waldecker et al. (2016, 2017), who predicted that phonons thermalize within a few picoseconds. The disagreement with their results is a consequence of the different theoretical description they proposed. They use a (non-microscopically) derived phonon branch dependence to account for the different strength in the electron-phonon interaction, and obtain the phonon-phonon coupling by fitting experimental data under the assumption of equal phonon-phonon interaction strength between branches. This assumption is however not justified as we have seen in Fig. 4, where the phonon lifetime is strongly mode dependent. Also, as we have seen in Fig. 2, the electron-phonon coupling constants do in fact strongly vary within the same phonon branch.

To provide the first ab initio description of laser-induced nonequilibrium dynamics in a metal, and to answer the question for how long the lattice is out-of-equilibrium, we have numerically solved our nonequilibrium rate Eqs. (2) and (3), considering phonon modes (enough to achieve good numerical convergence).

In Fig. 5 we plot the temporal evolution of the electronic and phonon mode temperatures for times up to 100 ps. The electronic temperature (blue line) increases rapidly, reaching up to 1575 K in 196 fs, followed by an exponential decrease with a decay time of about 1225 fs (light shaded area). The final electronic temperature reached (after 100 ps) is 366 K. We also show the phonon mode-dependent range of temperatures versus time (red area). The individual evolution of each phonon mode is left out for sake of simplicity. Initially, the maximum values in the range, which stem from optical phonon modes, increase, reaching a value of about 590 K at 1660 fs, followed by a slow monotonic decrease. On the other hand, the minimum values, which stem from acoustic modes close to the point, keep increasing very slowly, reaching a temperature of 313 K at 100 ps. Notably, the phonon mode temperatures cover a range of hundreds of Kelvins during several picoseconds, evidencing a nonequilibrium behavior. To estimate the weight of each phonon mode in the thermalization process and to avoid the singular behavior of some phonon modes, such as phonon modes close to , we also show the average phonon temperature (green line), calculated as the sum of phonon mode temperatures over . Although the average phonon temperature reaches an almost converged value very fast (within the first picoseconds), it is stunning to observe that this temperature still differs from the electronic temperature even at 100 ps, which confirms a continuing energy flow between the electronic and phononic systems. Moreover, the phonon modes temperatures cover an interval of about 50 K at 100 ps.

For comparison we have computed the electronic temperature evolution using the 2TM with the parameters obtained recently to simulate ultrafast demagnetization in FePt Mendil et al. (2014). The results, shown in the SM, are strikingly different. Not only does reach higher values (around 2500 K) in much shorter time (about 100 fs), also the final electronic temperature is higher, about 550 K. Additionally, the electron-lattice dynamics is much faster, reaching a common thermal equilibrium after only 1.5 ps. In contrast, our results evidence that the lattice remains out-of-equilibrium not only during short time scales after laser excitation but also on large timescales, and provide a clear example that a complete nonthermal modeling is needed to describe the out-of-equilibrium dynamics.

In conclusion, we have developed a nonequilibrium theory to describe the out-of-equilibrium dynamics triggered by ultrashort laser pulses. The model enables a fully ab initio determination of the out-of-equilibrium dynamics. Our simulations for FePt unambiguously reveal that, in contrast to previous understanding, the lattice is not in equilibrium even after 20 ps. As ultrafast nonequilibrium dynamics of materials is a strongly emerging research area and since our theory can provide a fully parameter-free description of the ensuing dynamics, we expect it to become a valuable tool for future modeling of ultrafast relaxation dynamics of laser-excited metals.

We thank H. Dürr and M. Bargheer for valuable discussions. This work has been funded through the Swedish research Council (VR), the K. and A. Wallenberg Foundation (Grant No. 2015.0060) and the Röntgen-Ångström Cluster. We also acknowledge support from the Swedish National Infrastructure for Computing (SNIC).

## References

- Fujimoto et al. (1984) J. G. Fujimoto, J. M. Liu, E. P. Ippen, and N. Bloembergen, Phys. Rev. Lett. 53, 1837 (1984).
- Beaurepaire et al. (1996) E. Beaurepaire, J.-C. Merle, A. Daunois, and J.-Y. Bigot, Phys. Rev. Lett. 76, 4250 (1996).
- Koopmans et al. (2010) B. Koopmans, G. Malinowski, L. F. Dalla, D. Steiauf, M. Fähnle, T. Roth, M. Cinchetti, and M. Aeschlimann, Nat. Mater. 9, 259 (2010).
- Liao et al. (2016) B. Liao, A. A. Maznev, K. A. Nelson, and G. Chen, Nat. Commun. 7, 13174 (2016).
- Hansteen et al. (2005) F. Hansteen, A. Kimel, A. Kirilyuk, and T. Rasing, Phys. Rev. Lett. 95, 047402 (2005).
- Ju et al. (1999) G. Ju, A. V. Nurmikko, R. F. C. Farrow, R. F. Marks, M. J. Carey, and B. A. Gurney, Phys. Rev. Lett. 82, 3705 (1999).
- van Kampen et al. (2002) M. van Kampen, C. Jozsa, J. T. Kohlhepp, P. LeClair, L. Lagae, W. J. M. de Jonge, and B. Koopmans, Phys. Rev. Lett. 88, 227201 (2002).
- Kim et al. (2012) J.-W. Kim, M. Vomir, and J.-Y. Bigot, Phys. Rev. Lett. 109, 166601 (2012).
- Kim et al. (2015) J.-W. Kim, M. Vomir, and J.-Y. Bigot, Sci. Rep. 5, 8511 (2015).
- Henighan et al. (2016) T. Henighan, M. Trigo, S. Bonetti, P. Granitzka, D. Higley, Z. Chen, M. P. Jiang, R. Kukreja, A. Gray, A. H. Reid, E. Jal, M. C. Hoffmann, M. Kozina, S. Song, M. Chollet, D. Zhu, P. F. Xu, J. Jeong, K. Carva, P. Maldonado, P. M. Oppeneer, M. G. Samant, S. S. P. Parkin, D. A. Reis, and H. A. Dürr, Phys. Rev. B 93, 220301 (2016).
- Rabitz et al. (2008) H. Rabitz, R. Vivie-Riedle, M. Motzkus, and K. Kompa, Science 288, 824 (2008).
- Harmand et al. (2013) M. Harmand, R. Coffee, M. R. Bionta, M. Chollet, D. French, D. Zhu, D. M. Fritz, H. T. Lemke, N. Medvedev, B. Ziaja, S. Toleikis, and M. Cammarata, Nat. Photon. 7, 215 (2013).
- Lindenberg et al. (2017) A. M. Lindenberg, S. L. Johnson, and D. A. Reiss, Annu. Rev. Mater. Res. 47, 426 (2017).
- Mitrano et al. (2016) M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S. Lupi, P. D. Pietro, D. Pontiroli, M. Ricco, S. R. Clark, D. Jaksch, and A. Cavalleri, Nature 530, 461 (2016).
- Kaganov et al. (1957) M. I. Kaganov, I. M. Lifshitz, and L. V. Tanatarov, Sov. Phys. JETP 4, 173 (1957).
- Anisimov et al. (1974) S. I. Anisimov, B. L. Kapeliovich, and T. L. Perel’man, Sov. Phys. JETP 39, 375 (1974).
- Allen (1987) P. B. Allen, Phys. Rev. Lett. 59, 1460 (1987).
- Stojchevska et al. (2014) L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D. Svetin, S. Brazovskii, and D. Mihailovic, Science 344, 177 (2014).
- Aoki et al. (2014) H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Rev. Mod. Phys. 86, 779 (2014).
- Carpene (2006) E. Carpene, Phys. Rev. B 74, 024301 (2006).
- Mueller and Rethfeld (2013) B. Y. Mueller and B. Rethfeld, Phys. Rev. B 87, 035139 (2013).
- Carva et al. (2013) K. Carva, M. Battiato, D. Legut, and P. M. Oppeneer, Phys. Rev. B 87, 184425 (2013).
- Baranov and Kabanov (2014) V. V. Baranov and V. V. Kabanov, Phys. Rev. B 89, 125102 (2014).
- Waldecker et al. (2016) L. Waldecker, R. Bertoni, R. Ernstorfer, and J. Vorberger, Phys. Rev. X 6, 021003 (2016).
- Waldecker et al. (2017) L. Waldecker, T. Vasileiadis, R. Bertoni, R. Ernstorfer, T. Zier, F. H. Valencia, M. E. Garcia, and E. S. Zijlstra, Phys. Rev. B 95, 054302 (2017).
- (26) The factor 2 stems from the spin degeneracy, but a more general form can be given.
- Ziman (1960) J. M. Ziman, Electrons and Phonons: The Theory of Transport Phenomena in Solids (Oxford University Press, Oxford, 1960).
- (28) Supplementary Material is available for this article.
- Gonze et al. (2009) X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken, F. Bottin, P. Boulanger, F. Bruneval, D. Caliste, R. Caracas, M. Cote, T. Deutsch, L. Genovese, P. Ghosez, M. Giantomassi, S. Goedecker, D. R. Hamann, P. Hermet, F. Jollet, G. Jomard, S. Leroux, M. Mancini, S. Mazevet, M. J. T. Oliveira, G. Onida, Y. Pouillon, T. Rangel, G.-M. Rignanese, D. Sangalli, R. Shaltaf, M. Torrent, M. J. Verstraete, G. Zerah, and J. W. Zwanziger, Comp. Phys. Commun. 180, 2582 (2009).
- Togo et al. (2015) A. Togo, L. Chaput, and I. Tanaka, Phys. Rev. B 91, 094306 (2015).
- Chase et al. (2016) T. Chase, M. Trigo, A. H. Reid, R. Li, T. Vecchione, X. Shen, S. Weathersby, R. Coffee, N. Hartmann, D. A. Reis, X. J. Wang, and H. A. Dürr, Appl. Phys. Lett. 108, 041909 (2016).
- Sun et al. (1993) C.-K. Sun, F. Vallée, L. Acioli, E. P. Ippen, and J. G. Fujimoto, Phys. Rev. B 48, 12365 (1993).
- Guo et al. (2001) C. Guo, G. Rodriguez, and A. J. Taylor, Phys. Rev. Lett. 86, 1638 (2001).
- Lisowski et al. (2004) M. Lisowski, P. A. Loukakos, U. Bovensiepen, J. Stähler, C. Gahl, and M. Wolf, Appl. Phys. A 78, 165 (2004).
- Lambert et al. (2014) C.-H. Lambert, S. Mangin, B. S. D. C. S. Varaprasad, Y. K. Takahashi, M. Hehn, M. Cinchetti, G. Malinowski, K. Hono, Y. Fainman, M. Aeschlimann, and E. E. Fullerton, Science 345, 1337 (2014).
- John et al. (2017) R. John, M. Berritta, D. Hinzke, C. Müller, T. Santos, H. Ulrichs, P. Nieves, J. Walowski, R. Mondal, O. Chubykalo-Fesenko, J. McCord, P. M. Oppeneer, U. Nowak, and M. Münzenberg, Sci. Rep. 7, 4114 (2017).
- (37) Our ab initio calculated lattice parameters of ferromagnetic FePt, Å and Å (), are in agreement with experimental values ( Å, Å, and ). The computed ground state magnetic moments , are also in good agreement with previous calculations and experiments (Oppeneer, 1998; Lyubina et al., 2006).
- Mendil et al. (2014) J. Mendil, P. Nieves, O. Chubykalo-Fesenko, J. Walowski, T. Santos, S. Pisana, and M. Münzenberg, Sci. Rep. 4, 3980 (2014).
- Lin et al. (2008) Z. Lin, L. V. Zhigilei, and V. Celli, Phys. Rev. B 77, 075133 (2008).
- Oppeneer (1998) P. M. Oppeneer, J. Magn. Magn. Mater. 188, 275 (1998).
- Lyubina et al. (2006) J. Lyubina, I. Opahle, M. Richter, O. Gutfleisch, K.-H. Müller, L. Schultz, and O. Isnard, Appl. Phys. Lett. 89, 032506 (2006).