Simulating electronically driven structural changes in silicon with two-temperature molecular dynamics

Radiation can drive the electrons in a material out of thermal equilibrium with the nuclei, producing hot, transient electronic states that modify the interatomic potential energy surface. We present a rigorous formulation of two-temperature molecular dynamics that can accommodate these electronic effects in the form of electronic-temperature-dependent force fields. Such a force field is presented for silicon, which has been constructed to reproduce the ab initio-derived thermodynamics of the diamond phase for electronic temperatures up to $2.5\mathrm{eV}$, as well as the structural dynamics observed experimentally under nonequilibrium conditions in the femtosecond regime. This includes nonthermal melting on a subpicosecond timescale to a liquidlike state for electronic temperatures above $\sim 1\mathrm{eV}$. The methods presented in this paper lay a rigorous foundation for the large-scale atomistic modeling of electronically driven structural dynamics with potential applications spanning the entire domain of radiation damage.

Figure 1

(a) The dependence of the $A$ and $B$ force field parameters on electronic temperature $T_e$. (b) The interatomic potential $U$ for a range of atomic volumes in the diamond phase, and for $T_e=0$, 1.5, and $2\mathrm{eV}$. The lines are the 2T-MD interatomic potential while the circles are the ab initio calculations. (c) The evolution of the nuclear temperature $T_n$ in response to a low-fluence excitation at $t=0$. The line is the 2T-MD calculation, the circles are experimental data points [1]. (d) The modified angle-dependence for $\lambda =0,0.2,0.4$, corresponding to a weakening of bond directionality. The inset shows the dependence of $\lambda$ on electronic temperature.

Figure 2

(a) Several normalized Bragg peaks decaying in response to a low-fluence laser pulse at $t=0$. The lines are predicted from 2T-MD while the circles are experimental measurements [1]. From top to bottom, the peaks are: (111), (220), (311), (331), (531), (620). (b) The normalized response of the (220) Bragg peak to a high-fluence laser pulse that nonthermally melts the crystal. The ab initio data (dashed line) was taken from [13] and the experimental data (circles) from Ref. [2].

Figure 3

The evolution of the centrosymmetry order parameter in response to an ultrafast laser pulse at $t=0$ for a range of absorbed fluences for simulations that (a) conserve energy and (b) do not conserve energy. A larger value corresponds to more order (centrosymmetry), but note the inverted axis.

Figure 4

The root-mean-square deviation of the atomic positions in response to a laser pulse at $t=0$ with the absorbed fluence labeled.

Figure 5

(a)–(d) The time evolution of the four Hamiltonian components in response to an ultrafast laser pulse at $t=0$ for a range of absorbed fluences $\mathrm{\Phi }$, as well as the (e) electronic temperature and (f) nuclear temperature. (g) and (h) are the same as (e) and (f), respectively, but for simulations that do not conserve energy.