ARTICLE pubs.acs.org/JPCC
Electron-Ion Coupling in Shocked Energetic Materials Evan J. Reed* Department of Materials Science and Engineering, Stanford University, Stanford, California 94304, United States ABSTRACT: The magnitude and role of electronic excitations in shocked energetic materials are studied theoretically using quantum molecular dynamics simulations. Focusing on the detonating primary explosive HN3 (hydrazoic acid), this work finds that the material transiently exhibits a high level of electronic excitation characterized by carrier densities in excess of 1021 cm3, or one electronic excitation for every eight molecules. Electronic excitations enhance the kinetics of chemical decomposition by ∼30%. The electronic heat capacity has a minor impact on the temperatures exhibited, on the order of 100 K. These simulations are performed using the self-consistent charge density functional tight-binding method (SCC-DFTB) combined with a new modification of a multiscale computational scheme for simulation of the coupling between electrons and ions in shocked matter.
’ INTRODUCTION When a compressive shock wave propagates through a material, high resulting temperatures and electronic structure changes can induce electronic excitations out of the ground state. It has been proposed that these electronic excitations play a role in a range of shock phenomena including the phase diagram and metallization of hydrogen,1,2 phase diagrams of water,3 helium,4 and carbon,5 and the detonation of energetic materials.68 In this work, the magnitude and role of electronic excitations are theoretically studied in a detonating chemically reactive explosive, liquid hydrazoic acid (HN3). HN3 is a highly sensitive liquid azide that is a chemically simpler analog of commonly utilized azides like lead azide and sodium azide. The calculations are accomplished by introducing an extension of a multiscale shock simulation technique (MSST)9 to enable treatment of electronion energy coupling in quantum molecular dynamics simulations of shock compression. The calculations show that electronic excitations can alter the potential energy surface upon which the classical ions evolve, leading to enhanced kinetics of the shockinduced decomposition. Also studied is the qualitative timeevolution of the electrical transport during the detonation. Relatively large levels of electronic excitation lead to a transition to an electrically conducting state for several picoseconds while chemical reactions occur. Energetic materials exhibit the phenomenon known as detonation where a net exothermic chemical reaction front propagates supersonically through the material as a shock wave. The temperatures and pressures in the shock are in the range of thousands of kelvins and tens of gigapascals, respectively, for condensed phase materials. The temperatures are high enough to generate a non-negligible level of thermal excitations of electrons if the bandgap of the material is sufficiently small. Both thermal and nonthermal mechanisms of electronic excitation have been proposed to play some role around shock fronts in detonation waves.68 In particular, it has been suggested that electronic excitations may be required to r 2012 American Chemical Society
achieve fast reaction time scales during detonation.6 Changes in electronic bandgap with pressure have been investigated in energetic materials using quantum simulation methods (e.g., ref 10), but these calculations neglect high-temperature bond distortions and chemical reactions. Molecular dynamics simulations of chemistry in shocked energetic materials have been performed with classical (analytical) interatomic potentials (e.g., ref 11), but such simulations do not account for electronic excitation. Quantum calculations of the first chemical steps of detonation in nitromethane using the MSST indicate the material exhibits a transient decrease in bandgap during chemical decomposition that could lead to a high level of electronic excitation.12 The direct method of simulating a shock wave propagating through a material involves propagating the shock down a long tube of atoms. Achieving shock propagation time scales >10 ps requires use of a sufficiently large number of atoms to preclude the use of quantum methods like density functional theory (DFT).1,13 Because quantum methods are required for treatment of electronic excitations under extreme conditions, an alternative approach must be utilized. The shock simulation method utilized here is an extension of the MSST.9,1416 Instead of the direct approach, the computational cell of MSST follows a Lagrangian material element through the shock wave, enabling simulation of the shock wave with significantly fewer atoms and enabling the use of quantum methods. MSST simulates steady shock waves (i.e., constant shock speed) by time-evolving equations of motion for the atoms and volume of the computational cell to constrain Special Issue: Chemistry and Materials Science at High Pressures Symposium Received: July 15, 2011 Revised: November 21, 2011 Published: January 03, 2012 2205
dx.doi.org/10.1021/jp206769c | J. Phys. Chem. C 2012, 116, 2205–2211
The Journal of Physical Chemistry C
ARTICLE
Figure 1. Schematic illustration of two equivalent versions of the electron-ion coupling scheme employed here. (a) Direct component of heat flow between electrons and a component that can be considered to occur via a thermal reservoir. The latter component is associated with changes in electronic single-particle states at constant temperature and occurs on the time scale of atomic vibrations. These two heat flow components are combined in scheme in (b), which is equivalent to panel a. In both cases, the heat exchange in the electron-ion system is adiabatic.
the shock propagation direction stress to the Rayleigh line and energy to the Hugoniot energy condition.9,14,15 For a specified shock speed, the latter two relations describe a steady planar shock wave within continuum theory. Uniaxial strain appropriate for a planar shock wave is imposed on the computational cell. The MSST has been demonstrated to reproduce accurately the sequence of thermodynamic states throughout the reaction zone of explosives with analytical equations of state, shock waves in amorphous Lennard-Jones, and amorphous Tersoff carbon.15,16 It has been used to simulate shocks in water,17 nitromethane,12 and other materials represented by DFT and tight-binding quantum energy models.18,19
’ ELECTRON-ION COUPLING Quantum simulations of materials at nonzero temperatures frequently use a thermostat that maintains the ion and electronic temperatures at a specified value of Te = Ti. The electronic temperature Te determines the population of single-particle electronic states fj(t) via the Fermi relation fj ¼ ðexp ðεj ðTe , tÞ μðTe , tÞÞ=kb Te Þ þ 1Þ1
ð1Þ
This assumption implies that the electrons are at equilibrium defined by a temperature Te. In general, the chemical potential μ(t), single-particle eigenvalues εj(t), and occupation numbers fj(t) are all time-dependent. In this work, we utilize a spinrestricted approach so there are two electrons per fully occupied state, fj = 1. Because eq 1 is a condition of constant temperature and not constant energy, energy is free to flow between the electronic degrees of freedom and an infinite reservoir at temperature Te. Energy flows into and out of the electronic degrees of freedom when the electronic levels are time-dependent during the MD simulation. The energy flow rate is given by Te(∂Se/∂t)|Te where the electronic entropy is Se ðfεi ðTe , tÞg, μðTe , tÞ, Te Þ ¼ kb
∑i ½fi ðtÞ lnðfi ðtÞÞ
þ ð1 fi ðtÞÞ lnð1 fi ðtÞÞ ð2Þ This approach is statistical in nature and aimed at obtaining ensemble averages by averaging over a sequence of configurations in time.
Complications arise in the case of a shocked material. Complications are primarily associated with the requirement to describe microscopic dynamics (e.g., chemical kinetics) rather than merely statistical information in an equilibrium state. In the shock case, the temperatures are time-dependent: Te(t) and Ti(t). Furthermore, the energy of the shocked material is not coupled to an infinite reservoir. The usual assumption of coupling to an energy reservoir becomes questionable on the short, picosecond time scales of shock-induced chemistry. Coupling to an energy reservoir introduces fictitious time scale(s) that may impact chemical kinetics of the shocked material. This section develops a scheme for the energy coupling of electrons and ions that occur in shocked materials. The goals are two-fold: (1) Allow for time-dependent temperatures Te(t) and Ti(t) and (2) employ an adiabatic treatment that does not couple the system to an external energy reservoir. In a shocked material, the electron and ion temperatures can initially be out of equilibrium behind the shock front; that is, Te(t) 6¼ Ti(t). Furthermore, the electrons and ions can each be in internal disequilibrium, where a temperature is ill-defined. In this work, the electrons and ions are each taken to be in sufficient internal equilibrium to have well-defined Te(t) and Ti(t). It is assumed that the electronic state occupations are given by the Fermi distribution at temperature Te(t), implying that the electrons are in equilibrium with each other at all times. Before discussing the general case where Te(t) 6¼ Ti(t), we will first consider the assumption Te(t) = Ti(t), valid on time scales longer than the electronphonon coupling time scale. This assumption is most valid during the multipicosecond and slower chemical transformations that occur after the initial shock compression. It may break down near the shock front where Ti(t) can potentially vary more rapidly than the equilibration time scale. The electron-ion disequilibrium case of Te(t) 6¼ Ti(t) can be incorporated at the expense of an additional empirical coupling parameter, to be discussed later. The following focuses on the case where Te(t) = Ti(t). The original MSST scheme exhibits a conserved energy quantity ~e0 given by 1 vs 2 ~v 2 1 p0 ð~v0 ~vÞ ð3Þ ~e0 ¼ ~e þ Q ~v_ 2 2 ~v0 2 s_ i |2 + (1/M)ϕ({rBi}), and dots denote with ~e = ∑i(1/2M)mi|A! time derivatives. Here Q is a mass-like parameter for the computational cell volume v, the total mass M = ∑imi is a sum 2206
dx.doi.org/10.1021/jp206769c |J. Phys. Chem. C 2012, 116, 2205–2211
The Journal of Physical Chemistry C
ARTICLE
of atom masses mi, and ϕ({rBi}) is the potential energy of the system, and we utilize the scaled coordinate scheme of ref 20 _s i where A is a matrix that contains r_ i t A ! where rBi AsBi and ! the computational cell lattice vectors and B s i is the scaled position of atom i. The constant shock speed is vs, and the initial computational cell volume and stress are v0 and p0, respectively. The shock speed is the constraint parameter that controls the amplitude of the shock. We wish to establish a new conserved energy quantity that incorporates the energy of the electronic degrees of freedom. When electronic excitations are considered, the potential energy ϕ({rBi},Te) becomes a function of the electronic temperature. There are two mechanisms of energy flow between the ions and electrons, illustrated in Figure 1. The first is the term discussed above, Te(∂Se)/(∂t)|Te , which is associated with changes in electronic entropy when single particle eigenstates vary in time. Heat flow can occur through this mechanism when the electronic temperature Te is constant. For example, the single-particle energy of a bond moves up and down as it undergoes vibration, resulting in a change of population and entropy of that state. A second mechanism of heat flow can occur when chemical evolution of heat changes Ti(t). Because it is assumed here that Te(t) = Ti(t), heat flows to electrons from ions at the rate dTe dTe ∂Se ðTe ðt Þ, t Þ Ce ¼ Te ∂Te dt dt t
Here Ce is the instantaneous heat capacity of the electrons. A modified MSST conserved quantity can be introduced that includes the two forms of heat flow between electrons and ions 0 0 1Z t 0 dTe ðt 0 Þ 0 ∂Se ðTe ðt Þ, t Þ 0 ~e0 dt Te ðt Þ jTe þ Ce ðt Þ M t0 ∂t 0 dt 0 ð4Þ
times. The remainder of this work focuses on the case of Te(t) = Ti(t) because the disequilibrium case comes at the expense of an additional empirical parameter α, which is not known and is furthermore likely to depend on temperature, pressure, and chemical composition. In the absence of accurate knowledge of this parameter for HN3, we make the assumption of instantaneous equilibrium. This approximation implies the assumption that the time scale of the chemistry is longer than the electronphonon coupling time scale under these conditions. The MSST equations of motion of the atoms must be altered to incorporate heat flow between electrons and ions. The original equations of motion for the computational cell volume and atoms are :: Q ~v ¼
ð7Þ
! s_ i jA! s_ j2
ð8Þ
dA
dϕ
2 :: 1 1 ∂ϕ μ~v_ 1 _ ! ! _ A G G si þ M si ¼ mi ∂B ri mi~v
∑j
j
where μ is an empirical macroscopic viscosity term and G ATA.16 Uniaxial strain of the computational cell is utilized because the shock is assumed to be planar and (dA)/(dv) is a matrix corresponding to the direction of the uniaxial strain; A is allowed to change with only one degree of freedom. These equations of motion can be derived from the Lagrangian form of eq 3 when μ = 0. The last terms containing μ in these equations couple kinetic energy of the volume degree of freedom to the individual atoms, acting as artificial viscosity beyond that intrinsically provided by ϕ({rBi}). To account for heat exchange with electrons, an extra term is included in the equation of motion of the atoms
Because Te ðt 0 Þ
~v_
v2
∑i miA!s_i 3 dv !s_i dv ~vs20 ð~v0 ~vÞ p0 μ~v
2 ! :: 1 1 ∂ϕ Mμ~v_ dSe s_ i ! si ¼ Te s_ i þ A G1 G_ ! 2 ! _ mi ∂B ri ~ v dt mi jA s j j
0 ∂Se ðTe ðt 0 Þ, t 0 Þ dTe ðt 0 Þ 0 0 dTe ðt Þ jTe þ C ðt Þ ¼ T ðt Þ e e ∂t dt 0 dt
∑j
ð9Þ
the conserved quantity can be written in a more familiar form ~e0
1Z t 0 dSe ðt 0 Þ dt Te ðt 0 Þ M t0 dt 0
ð5Þ
Note that if Te(t) is constant in time, then the conserved quantity reduces to the usual free energy ~e0 Te~Se(t), as might be expected. In general, however, Te(t) is time-dependent, and expression 5 is the relevant conserved quantity. So far, instantaneous equilibrium between electrons and ions has been assumed, Te(t) = Ti(t). A related coupling scheme can be described for the case of Te(t) 6¼ Ti(t), where the ions and electrons are not assumed to be in equilibrium. A coupling parameter α can be defined by expressing the heat flow as Te(t)(dSe(t))/(dt) α(Ti(t) Te(t))Ce(t). The electronic temperature changes according to Te ðtÞ ¼
αTi ðtÞCe ðTe , tÞ dSe ðTe , tÞ þ αCe ðTe , tÞ dt
ð6Þ
which reduces to the case of Te = Ti in the large α limit. The isentropic case of (dSe/dt)(Te,t) = 0 is established in the limit of small α, where Te is adjusted to set (dSe)/(dt)(Te,t) = 0 at all
!
where Te is taken to be the instantaneous ion temperature Te ðtÞ ¼ Ti ðtÞ
1 3NkB
∑j mjjA!s_j j2
ð10Þ
The volume equation of motion remains unchanged. The addition of the extra term in eq 9 provides::an extra force on the ions B F j that v j = ∑jmjA! s_ j = Te(dSe/dt). ThereF j 3B s j 3 A! yields a power ∑jB fore, the energy condition eq 5 is conserved; that is, changes in ϕ({rBi},Te) associated with electronic excitation are compensated by the entropy term in eq 9. Equations 7 and 9 can be integrated using methods for the Langevin equation for which iterative Verlet-based time-reversible methods have been described.21 By relaxing overall time reversibility, a noniterative Verlet-like method can be used in conjunction with expressions for evaluating the quantities ! s_ i (t), _~v (t), and (dSe/dt)(t) at time t.22 A similar scheme has been used for the calculations in this work. At the expense of sacrificing time reversibility, the relation dSe ðtÞ ¼ 3Se ðtÞ 4Se ðt ΔtÞ þ Se ðt 2ΔtÞ þ OðΔt 2 Þ dt 2207
ð11Þ
dx.doi.org/10.1021/jp206769c |J. Phys. Chem. C 2012, 116, 2205–2211
The Journal of Physical Chemistry C
Figure 2. Time dependence of temperature and conserved energy quantity ~ε0/3NkB for a 7 km/s shock in HN3. Here the conserved energy quantity is normalized by 3NkB, a rough approximation for ionic heat capacity, to give an approximate indication of temperature drift associated with integration errors. The inset shows the Rayleigh line target stress for this shock speed along with the simulation trajectories for the intrinsic stress and stress including empirical term (μ). These plots demonstrate adherence to the Hugoniot and Rayleigh constraints. See the text for details.
ARTICLE
Figure 3. Calculated temperature, stress, and volume versus time behind shocks propagating at 7.5 and 6 km/s (near detonation speed) through HN3.
is acceptable for evaluating the electronic entropy time derivative at time t in eq 9. Additional details regarding MSST equations of motion are in ref 16 and references therein.
’ APPLICATION TO LIQUID HN3 As an application of this method, this section considers the explosive liquid azide HN3 compressed by a shock propagating near the natural detonation speed of the material. The simulations here utilize the self-consistent charge density functional tight binding (SCC-DFTB) method in conjunction with the MSST.23,24 The Hamiltonian parameters utilized are designed for periodic boundary conditions (pbc-0-1 and pbc-0-3) and are available at http://www.dftb.org/parameters/download/. Electronic properties predicted by DFTB have been found to provide reasonable agreement with DFT calculations of the energetic material nitromethane under pressure.10,25 The first chemical reactions in the decomposition of nitromethane near detonation conditions have been found to be identical for DFTB and DFT.12,26 The energy of formation of HN3 (without zero-point energy) using the DFTB scheme is 0.083 H, slightly less than the value of 0.092 H calculated using DFT at the B3LYP/ 6-31G* level and 0.114 calculated at the QCISD/cc-PVTZ level (quadratic configuration interaction with single and double substitutions, ref 27). Whereas DFTB is an approximate method, it has the potential to provide insight into size and time scales beyond DFT and provide a description of the electronic states. Simulations in this work utilized an orthorhombic computational cell containing 64 molecules and employed periodic boundary conditions with the gamma point utilized for Brillouin zone integration. The molecular liquid was equilibrated with computational cell lattice vectors a = b = 12.79 Å and c = 25.58 Å, where uniaxial compression occurs along the c axis. The latter give an initial density equal to the experimentally observed density of 1.092 g/cm3. The MSST equations of motion utilized computational cell mass-like parameter Q = 1012 kg2/m4 and explicit macroscopic viscosity μ = 102 kg/m/s.
Figure 4. Density of states at different times during the 6 km/s shock. The Fermi energy is at zero. One picosecond averages around 0, 5, and 30 ps are shown, in addition to an average over the entire simulation (red curve). The DOS at the Fermi energy exhibits a transient increase during peak chemical reactivity (blue curve).
Figure 2 shows how the Hugoniot and Rayleigh constraints are satisfied for a 7 km/s shock speed simulation in HN3. An integration time step of 0.25 fs with a SCC convergence parameter of 104 was utilized with the dftb+ code to generate this plot. Plotted is the conserved energy quantity ~ε0/3NkB given by the right-hand side of eq 3 with the equations of motions for volume and atoms taken to be eqs 7 and 9, respectively. Here the conserved energy quantity is normalized by 3NkB, a rough approximation for ionic heat capacity, to give an approximate indication of temperature drift associated with integration errors. The inset shows the trajectories of the intrinsic stress (first two terms of eq 7) and the intrinsic stress plus empirical viscous component (including last term of eq 7) overlaid on the Rayleigh line (remaining terms on right-hand side of eq 7). These trajectories show that the empirical viscosity plays a role during the initial compression but plays minimal role during the subsequent chemistry and volume expansion. The simulations start 2208
dx.doi.org/10.1021/jp206769c |J. Phys. Chem. C 2012, 116, 2205–2211
The Journal of Physical Chemistry C
Figure 5. Time dependence of the hole and electron densities during a 6 km/s shock. Holes dominate during the time of peak chemical reactivity because the Fermi level is located closer to the valence states than conduction states. (See Figure 4.) At their peak, these carrier densities correspond to ∼0.1 carriers per molecule.
Figure 6. Time dependence of TeSe normalized by 3NkB, a rough approximation for ionic heat capacity. This temperature is an approximate measure of the temperature shift of the shocked material due to electronic heat capacity. For both 6 and 7.5 km/s shocks, this parameter exhibits a peak when the density of states around the Fermi energy is maximal, followed by a smaller value at later times. The magnitude of these temperatures is a few percent of the 20006000 K temperatures exhibited by these materials. (See the temperatures in Figure 3.)
at V/V0 = 1 and compress to the left, then expand more slowly to the right along the Rayleigh line as chemistry occurs. Figure 3 shows temperature and stress versus time behind the shock front for a 6 and 7.5 km/s shock propagating through HN3, near the 7.1 to 7.6 km/s range of experimentally measured detonation speeds.28,29 The volume decreases rapidly during the initial compression, followed by a slower volume increase as chemistry occurs. The slower expansion is accompanied by a temperature increase as heat is evolved from the reactions. Figure 4 shows the density of states (DOSs) at different times during the 6 km/s shock. The material initially has an electronic gap >2 eV, but the gap has disappeared after 5 ps during the peak chemical reaction period (blue curve). At late times (30 ps, green curve), the DOS at the Fermi energy decreases again and a gap opens up. Similar nonmonotonic temporal time evolution of the DOS near the Fermi energy was also observed in nitromethane under detonation conditions in previous work.12 The high DOSs at the Fermi energy around 5 ps combined with multithousand degree temperatures generate electronic excitations. Five shows the concentration of holes and electrons in the material. The hole density is taken to be Fh = (2/v(t))∑i
ARTICLE
Figure 7. Time dependence of dominant species populations during 6 km/s shock, near detonation speed. Thin and thick lines correspond to the case of 200 K fixed electronic temperature (Te) and the case where electron (Te) and ion (Ti) temperatures are set equal. Electronic excitations lead to enhanced kinetics in the case of Te = Ti.
(1 fi(t))Θ(2fi(t) 1), and electron density is taken to be Fe = (2/v(t))∑ifi(t)Θ(1 2fi(t)), where Θ is the Heaviside function and the sum is over all single-particle electronic states i. Holes dominate during the time of peak chemical reactivity because the Fermi level is located closer to the valence states than conduction states. (See Figure 4.) During peak times, these carrier densities correspond to ∼10 holes in the computational cell, or ∼1 hole for every 8 molecules. To the author’s knowledge, the carrier concentrations observed during detonation have not been experimentally measured. An experimental observation has been reported of a positively charged region propagating in front of the detonation wave in nitromethane.30 It is not inconceivable that such a region may be created by the diffusion of holes in front of the shock. Figure 6 shows time dependence of the temperature Te(t)Se(t)/ 3Nkb which is a rough measure of the temperature shift of the system due to the heat capacity of the electrons, taking an approximate ionic heat capacity of 3Nkb. For shocks at 6 and 7.5 km/s speeds, this parameter exhibits a peak when the DOSs around the Fermi energy is maximal followed by a smaller value at later times. The magnitude of these temperatures is a few percent of the 20006000 K temperatures of these materials, shown in Figure 3. This indicates that electronic heat capacity alters temperatures by a few percent along the decomposition pathway. Electronic heat capacity here is on the order of other neglected effects like the quantum nature of the nuclear vibrations that can shift shock temperatures by a few hundred degrees kelvin.31 Figure 7 shows a comparison of the time dependence for dominant molecular populations for calculations with and without the electron-ion coupling scheme. Molecules were defined by interatomic bonds
