{'fraction_non_alphanumeric': 0.025953639265344566, 'fraction_numerical': 0.006854691015645439, 'mean_word_length': 4.651094027202839, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Molecular Dynamics Simulations of Detonation Instability

After making modifications to the Reactive Empirical Bond Order potential for Molecular Dynamics (MD) of Brenner et al. in order to make the model behave in a more conventional manner, we discover that the new model exhibits detonation instability, a first for MD. The instability is analyzed in terms of the accepted theory.


I. INTRODUCTION
There are obvious problems with the 1D theory of detonation. Foremost is that it is impossible to create truly 1D detonation experiments. In many cases experiments come close, and their results can be extrapolated to match theory. But for certain classes of explosives, namely gas phase ones, the shape of the shock front, no matter how one-dimensional the experimental setup is, is intrinsically 3D with ripples. These detonations leave patterns in the soot on the walls of shock tubes, indicating the incidence of a section of the front with high vorticity with the wall. Such detonations are said to be unstable. See for example Chapter 7 of Fickett and Davis's book [1] or Austin's dissertation [2].
The first truly successful work in theoretically analyzing this phenomenon was done by Erpenbeck in a series of papers published throughout the sixties [3,4,5,6] (summarized in chapter 6 of ref. [1]) in which he describes a Laplace transform procedure for determining stability for a given parameterization of the equation of state (EOS) and reaction rate of a material. In 1964 he applied his formulation to 1D instability for a variety of parameterizations [4], and in 1966 he expanded to transverse perturbations (multi-D) [5]. In 1970 he summarized his decade of work [6].
Although Erpenbeck's method was the only one at the time to definitively show whether a parameterization was stable, it was not widely used for reasons mentioned by Lee and Stewart [7]. The stability of a detonation based on a given parameterization had to be figured point-wise and the stability boundary interpolated. Lee and Stewart introduced a normal mode analysis via shooting method [7] that is easier and faster than Erpenbeck's Laplace transform method. The normal mode analysis has fewer drawbacks about in which regimes it is useful, namely underdriven detonations.
After the introduction of normal mode analysis, there was an expansion in detonation instability research. Notable is the work of Bourlioux and Majda [8,9], who were the first to compare the results of the normal mode analysis with hydrodynamic detonation simulations in 2D. Short and Stewart extended the work of Lee and Stewart to transverse perturbations [10].
Most analysis and simulation have been on square detonations (instantaneous shock rise, finite induction zone, and instantaneous heat release, typical of detonations of high activation energy [1]) with a polytropic gas EOS and an Arrhenius reaction rate with one-step reactions (although other rates, numbers of reaction steps, and different EOSs have been ex-amined [11]). The main parameters of concern in these analyses are the activation energy, exothermicity, adiabatic gamma, and degree of overdrive. Short gives a good summary of the research done since normal mode analysis in the area of detonation instability [11]. Two notable observations are that all of the simulations are done by continuum hydrodynamics codes and that condensed phase high explosives are, for the most part, stable. Instability has never been shown in the current range of limited molecular dynamics (MD) studies, which make no assumptions about the EOS and reaction rate.
In the preceding paper [12] we justify modifications to the functional form of a popular MD model (ModelI) for high explosives (HEs) (see Eq. 1 and Table I in the errata of Brenner et al. [13]). The result is a new potential (ModelIV), which behaves in a more conventional manner. We notice, however, that the measured shock velocity is noticeably faster than the Chapman-Jouguet (CJ) value for an underdriven detonation. Where the preceding paper [12] focuses mainly on the difference in thermo-chemical properties of the two models, this current work is concerned with the results of non-equilibrium MD (NEMD) simulations of the ModelIV potential, which will demystify the shock velocity discrepancy. In Sec. II we explain our choice of dissociation energies. In Sec. III the detonation instability is measured and analyzed in terms of theory.

II. PARAMETERIZATION
Two types of 2D simulations using the ModelIV model are performed. The preceding paper [12] discusses the results of microcanonical simulations for seeking the CJ state and conducting Arrhenius cookoffs. There a parameterization of ModelIV in which the metastable dissociation energy D AB e = 1.0 eV and the stable dissociation energy D AB e = 5.0 eV such that the exothermicity Q = D AA/BB e − D AB e = 4.0 eV is used, whereas with ModelI both D AB e and Q were never varied away from the REBO defaults [13] at once. This parameterization was initially chosen to create a self-propagating detonation front. In this section the results of NEMD simulations are discussed, in which supported and overdriven detonations are modeled. The reason for the choice of parameterization is discussed.
Initially, the same bonding energies as with ModelI, D AB e = 2.0 eV and D AA/BB e = 5.0 eV, are used. The first characteristic discovered is that the new material is insensitive to shock initiation. In fact, initiation requires a rigid piston driving at a velocity of about 4.5 km/s, an overdriven state. Since it is presumed that this new material has a high activation energy, but sufficient exothermicity to maintain detonation once initiated even if unsupported, the piston's velocity is backed off to zero. The result is that the detonation quickly quenches once the rarefaction wave that emanates from the decelerated piston catches up to the reaction zone. It seems, then, that we do not have a true HE.
It should be noted that even with a strong overdriving piston, the material at the front demonstrates a resilience not found in ModelI. Where ModelI succumbs to nearly instantaneous dissociation upon compression by an unsupported detonation, ModelIV demonstrates uniaxial compression in a finite induction zone, characterized by vertically aligned (detonation travels in the horizontal direction) reactant dimers before dissociation.
Since it is assumed that the exothermicity is sufficient, the parameters for a second simulation are next set to D AB e = 1.0 eV and D AA/BB e = 4.0 eV such that Q still is 3.0 eV. This detonation, which, as with the previous parameterization, is initiated by an overdriving piston that was backed off to zero, eventually fails too (see pane a in Fig. 1). Upon inspection of Fig. 2, one will notice that the front seems to have a standing wave shape, but Fig. 3 shows that there are symmetric traveling waves, probably emanating from either side of the same event.
Structure in the front is a tell-tale sign of instability in the detonation and violates the 1D assumption of the Zel'dovich, von Neumann, and Döring theory [1]. Such structure has been theorized [3] and shown through experiments and hydrodynamic simulations-mostly in reactive gases- [1] to be caused by transverse waves crossing in the reaction zone. The keystone shape in the λ contour in the first row of Fig. 2 is also typical of such unstable detonations in gases [2].
The number density (η) contours in Fig. 2 are plotted in order to see these postulated transverse waves, but the results are inconclusive. Notice that the degree of reaction (λ) and η seem to have an inverse relationship within the reaction zone. Behind a trough in the shock front, η is high and λ is low. Eventually the material there reacts. The resulting energy release pushes the products out of the area thus surging the trough forward.
Since it is hypothesized that the distance between the periodic boundaries in the transverse direction dictates the wavelength of what seemed to be a standing wave, a third simulation with all of the same parameters, but with a sample that is twice as wide (192 lattice cells as opposed to 96), is started. This simulation quenches quickly (see pane b in Fig. 1) so a fourth is performed with an initiating piston that starts out moving fifty percent faster, but that, too, does not sustain the detonation. If there are transverse waves, they may help to sustain the detonation in the thinner simulation. Increasing the distance between the boundaries may have dropped the frequency of crossings of the transverse waves below a critical limit. This also implies that the main shock front is insufficient to sustain detonation.
In order to make the detonation be self sustaining, a third parameterization in which D AB e = 1.0 eV and D AA/BB e = 5.0 eV, such that Q = 4.0 eV, is tried. So far this simulation, which, as with the preceding, is initiated with an overdriving piston that backs off to zero, has not shown signs of failure (see pane e in Fig. 1). It is with this parameterization that this research proceeds. This simulation is as thick as the thicker of the two previous (192 lattice cells thick). In this simulation the front has an asymmetric structure (see Fig. 4).

III. DETONATION INSTABILITY
Notice the distinct cellular structure in the λ contours, especially in the latter frames of Fig. 5. This is typical of unstable detonation and is seen in soot prints left on the walls of rectangular tubes used in detonation experiments and seen in certain properties in hydrodynamic simulations of the reactive Euler equations [1]. Where experiments using condensed phase explosives show a break down of this regular cell structure, ModelIV has a clear cell structure. There are many differences between ModelIV and real condensed phase explosives, particularly in the current implementation, which is performed on perfect crystals while real materials have a high concentration of defect structures. Any of said differences could mean the difference between cellular structure and random structure, but ModelIV suggests there is more of a reason for the real break down than merely being in the condensed phase.
From the η contours (Fig. 5), one can clearly see transverse waves (perhaps so here and not in Fig. 2 because here the resolution is finer and Q is greater so that the gradients at the fronts of the transverse waves are also greater). The sequence of events displayed in Fig. 5 can be matched to the accepted theory, of which Fickett and Davis give a good summary in Chapter 7 of their book [1]. Following their example, start with a flat shock wave, followed by an induction zone and then a reaction zone, a 1D situation. If a point "microdetonation" occurs within the induction zone, it will produce a shock wave of its own, the shape of which will be skewed by the rarefaction and flow already present. It will have an arched shape. Since the shock stability conditions state so, the forward facing part of the microdetonation's shock will overtake the original flat front causing it to bulge forward.
If there are two microdetonations, the transverse facing parts of their shock fronts will eventually meet. That intersection will move forward toward the shock front. Once the angle between the colliding transverse waves increases sufficiently, a Mach stem is formed. If conditions are right, this occurs within unreacted material, and their collision causes reaction, which surges the Mach stem forward of the incident shock front and intensifies the transverse waves as they pass through each other.
As the Mach stem ages, the induction zone behind it lengthens because of rarefaction. The sections of transverse waves near the shock front move as detonation fronts through the rest of the induction zone, consuming post-shocked yet unburnt material on either side of the Mach stem. The transverse waves are not sustainable without fuel. If they again collide behind an old Mach stem (or any place with a long induction zone), they will have plenty of fuel to cause the formation of another rapid burn and forward surge. If a transverse wave encounters a newly formed Mach stem, behind which there is very little unburnt fuel, it will weaken. This gives spacing criteria for surviving transverse waves: the distance between them must be great enough that sufficient growth of the induction zone can take place between passings and small enough that new waves do not have sufficient space in which to grow and consume the fuel between existing waves. Since the NEMD simulations use periodic boundary conditions, the wave spacing is at most the simulation thickness. Since it was shown that these unsupported detonations cannot sustain themselves without the transverse waves, their ability to propagate depends on the simulation width as is evident in Fig. 1.
The intersection of a transverse wave with the shock front is called the triple point. It is the passing of a triple point that etches cell patterns into the soot in shock tubes. It is also the source of the cell patterns seen in our λ contours. Experiments cannot see cell patterns in chemistry, at least not in OH fluorescence from detonating O 2 + H 2 [2]. Since, in thriving transverse waves, the triple point is crossing unreacted material and causing it to detonate, the shock front behind the triple point surges forward. The track of a triple point is thus seen as a traveling trough in contour plots of the relative shock position.
If a transverse wave encountered only constant state material, it would reach a steady state. But, since there are variations by collisions with other transverse waves, passing through its own wake via reflection (or, in the case of the NEMD simulation, periodic boundaries), Mach stem formation, and thermal fluctuations, etc. there can be a chaotic and stochastic quality to the behavior and number of transverse waves, the amount of which should depend on the transverse waves' sensitivity to such variations as well as the dimensions of the simulation. This dependence is evident in the degree of regularity in cell patters and and the evolution of the shock surface (see Fig. 4).
The example above is what one sees in Fig. 5 (and in the supplemental videos [15]). In the first η frame (f) one can clearly see two transverse waves receding after the formation of a Mach stem, which has surged forward. The triple points of the transverse waves are headed toward unreacted material. Notice in the analogous λ contour (a) that the induction zone is short behind the Mach stem. A final observation for this frame is that there is a plume of reactant that escapes detonation. It burns up more slowly as it flows downstream of the shock front. This is typical of unstable detonations [2].
In the second row of Fig. 5 (frames b and g), the triple points of the two transverse waves are about to collide in the lower half, where material is unreacted. The Mach stem is beginning to rarefy. We trace out the relevant lines from this row in Fig. 6. In this Figure one can make out the two triple In row three of Fig. 5 (frames c and h), the material has combusted after the crossing of the triple points, and the front has surged forward, a state seen in the first row. A second Mach stem has formed while the original one has further rar- efied, providing fuel for the yet again strengthened transverse waves.
In rows four (d and i) and five (e and j) another Mach stem is formed, but the second one has not rarefied much yet. Because of a lack of fuel, the downward passing wave eventually disperses, as is evident in Fig. 4 by the disappearance of the dark line of negative slope shortly after 50 ps.

IV. CONCLUSION
For ModelIV the activation energy is higher than with Mod-elI. The fact that instability appears in ModelIV but not Mod-elI follows the trend that stability tends to increase with decreasing activation energy [10]. As far as we know, this is the first time such instability has been demonstrated using an MD model.
The hope is that the ability for MD to model instability will eventually reveal much about the phenomenon. It can be helpful because it needs no reaction rate or equation of state as input parameters and thus cuts down on some of the parameter space to be explored. For the former reason, it can be especially useful for the study of condensed phase high explosives since it is difficult not only to know their multiphase EOS throughout the reaction zone, but also to make real measurements in the reaction zone. Although this may not be too important for ZND like detonations, in which the unsupported shock velocity and the CJ state depend on the initial state and that of the products, it does become relevant in weak [1] detonations, in which the sonic point is at a state of incomplete reaction and which are a starting basis for instability analysis.
If one thoroughly probes the EOS and the reaction rate of ModelIV, normal mode analysis could confirm its predictions. The scaling of Short and Stewart [10] or one based on the CJ state may allow for comparison to real explosives.
It should be noted that the cellular structure found with ModelIV is not typical of most condensed phase explosives in either experiment, continuum modeling, or normal mode analysis, yet this is not disheartening because of the aforementioned problems each of these techniques has with such explosives. Experiments may not be able to resolve them. Hydrodynamic codes and normal mode analysis may be using the wrong EOS or reaction rate, or ModelIV may be missing some characteristic of condensed phase explosives that tends to stabilize the detonations thereof [11].
Many of these cited studies go into much detail about the parametric regimes in which their explosives are stable or not. They vary the activation energy, Q, the amount of overdrive of a piston, or the frequency or wavenumber of a perturbation for instance. Such a study would be worthwhile using ModelIV, but it is beyond the scope of the current research. However, this work has raised this as a significant issue, and shown that it can be studied with this approach.
Of greatest significance is that this work has shown condensed phase systems with realistic molecular properties can sustain detonations with 2D instabilities. The net result is that the systems propagate with velocities somewhat higher than CJ and that their reaction zone structure is convoluted with the structure of the transverse waves. All analyses of condensed phase HEs presume that they are 2D stable, as this is the only manner in which a complete analysis can be achieved. An example is the Cheetah EOS [14] which presumes that the observed detonation velocity is in fact the true CJ velocity. The current work raises the possibility that these values may be several percent too high. Other features such as corner turning depend on the width of the effective reaction zone, where one must now consider how transverse waves may effect this. Consequently, it is important to have developed a tool with which to quantitatively explore these features.