High-pressure deformation and amorphization in boron carbide

Hard materials are often classified based on their load-independent Vicker’s hardness.¹ While some of them are naturally (intrinsic) hard, others can be processed to possess high hardness (extrinsic). Materials having diamondlike structure (e.g., natural diamond, cubic boron nitride, etc.), called “ultrahard,” are known to have the highest intrinsic hardness in the range of 80–200 GPa, while some synthesized allotropes of carbon, e.g., fullerite $C60$⁠, possess extrinsic hardness of up to 300 GPa.² These compounds consist of strong covalent bonds and a dense network of carbon atoms, or its neighbors in the periodic table. Diamondlike structures, as well as fullerite $C60$⁠, have low resistance to heat, as their structures collapse at temperatures of about $1000°C$⁠. They are sensitive to chemicals and cannot be mass-produced. Next in hardness are materials possessing a boron-rich icosahedral structure, such as boron carbide, boron suboxide, different allotropic forms of boron, and recently developed boron-aluminum-magnesium (BAM) materials, which as a group, possess hardness in the range of 25–45 GPa.^3–6 With lower density (⁠$∼2.5g/cm2$⁠) than diamondlike materials (⁠$∼3.5g/cm2$⁠), high thermal and chemical stability, and ability to be mass-produced, along with a combination of favorable thermomechanical properties, they are suitable for many mechanical, electrical, electronic, magnetic, thermal, nuclear, and optical applications.⁷ A common drawback in structural applications of most boron-rich icosahedral solids is “pressure-induced amorphization,” propensity for catastrophic structural collapse at an atomic level, at conditions of high pressures. This susceptibility, which limits their applicability to resist high loads,⁸ has been the subject of extensive investigation for over two decades.⁹ The common underlying drive in pursing fundamental studies of boron carbide as an ultrahard material rests upon understanding the mechanisms of amorphization, as well as the manifestations thereof.

The basic structure of boron carbide [space group 166, (⁠$R3¯m$⁠)] contains a 15-atom unit cell having a boron-rich icosahedron and a 3-atom intericosahedral chain. Boron carbide is known to possess different polytypes and polymorphs based on varying carbon content and site-specific placement of carbon atoms^10,11 in the icosahedra. Due to the presence of closely resembling structures, which causes local changes in stoichiometry, researchers argue the overall stoichiometry of boron carbide to be $B4.3C$⁠.¹² Kunka et al.¹³ devised a spectral strategy to identify the presence of different polymorphs to quantify their volume fractions in a fabricated sample of boron carbide. Given the strong dependence of properties on polymorphic variations, a new nomenclature has been developed to effectively include a site-specific description.¹¹ Upon application of hydrostatic loads, the 3-atom chains in boron carbide are susceptible to bending, whereby the central atom gets pushed out of the collinear arrangement, which, according to some researchers, represents the onset of structural instability.¹⁴ This interstitial space, where the central atom pushes into, has been named “cage space.”⁹ 

At high stresses beyond the Hugoniot elastic limit (HEL), boron carbide exhibits structural collapse through pressure-induced amorphization, which is observed in shock tests,¹⁵ indentations,^16,17 and diamond-anvil cell experiments.¹⁸ Evidence for amorphization is found in TEM images,^15,17 wherein amorphous islands have been observed, spanning 1–2 nm in width and over a length of 100 nm,⁹ as shown in Fig. 1(a). The most interesting observation lies in the Raman spectra of its amorphized state. While unstressed boron carbide shows spectral peaks below $1200cm−1$⁠, pressure-induced amorphized boron carbide is known to exhibit new spectral features, including broad peaks situated near $1330cm−1$ (or $1340cm−1$ according to Refs. 19 and 20 and so referred to in the remainder of the paper), $1520cm−1$⁠, and $1810cm−1$⁠;^16,21,22 see Fig. 1(b). In addition to these dominant characteristics in the amorphous spectrum, we also observe a small right shift of the main crystalline peak²² [see Fig. 1(b)]. There have been a few attempts aiming to resolve the appearance of new peaks in Raman spectra.^16,19 These efforts claim the presence of new products of amorphization including amorphous carbon, graphite, $α$-boron (⁠$B12$⁠), and (⁠$B12$⁠)CCC,¹⁰ typically arising from bending of chains in boron carbide.²³ However, due to the absence of alternate experimental characterization tools at the nanoscale, existing theories bear the burden of being mere speculations. To date, there is no resolution on why boron carbide exhibits new spectral features of a special kind, upon amorphization. In order to bring out the true potential inherent in boron carbide to exist as the best candidate for ultrahard materials, it is important to discern the right manifestation of amorphization. Since this critical phenomenon is sudden and occurs at high pressures, there is a lack of thorough understanding of mechanics of material behavior, especially because the driving phenomena originate at the atomic scale.

The goal of contemporary research is to understand and address underlying mechanisms of amorphization in icosahedral boron-rich materials and ultimately design microstructures to improve their hardness close to the range of diamondlike structures. In the present article, we develop a new rationale of amorphization in boron carbide, identifying links between thermodynamic conditions of amorphization and the unique Raman spectrum exhibited by pressure-induced amorphized boron carbide. We perform density functional theory (DFT) simulations to obtain Raman spectra for boron carbide over a range of high pressures. Next, results from molecular dynamics (MD) simulations of volumetric compression are presented, which use a newly developed ReaxFF potential for boron carbide. These results are used to establish a thermodynamic framework of amorphization, identifying the key impact of pressure and temperature. This is followed by the derivation of analytical estimates of atomic-scale residual pressure exerted by amorphous islands. Finally, by proposing a mechanistic layout of residual stress, we demonstrate how DFT-evaluated Raman spectra, in the stress range 60–80 GPa, fit most of the new spectral features shown by the experimentally obtained spectrum of amorphized boron carbide. For completeness, the article also includes critiques of various theories of pressure-induced amorphization prevalent in the literature pertaining to icosahedral boron-rich ceramics, along with their effectiveness in unifying all experimental observations to date. Results presented in this paper also indicate the importance of multiscale modeling, especially the connection between the origin of amorphization at the atomic scale and its manifestation in Raman spectra in the micrometer scale.

In this section, we present a brief history of attempts to understand amorphization in boron carbide through experimental and computational approaches.

Discovery of the amorphization phenomenon in boron carbide dates two decades back, when Grady²⁴ performed shock-impact experiments on boron carbide and observed a catastrophic loss in strength above the Hugoniot elastic limit (HEL), concomitant with near-fluidlike behavior. He coined the term “postyield stress softening” to differentiate boron carbide from silicon carbide, which rather exhibits near-metal-like hardening behavior, typical of other structural ceramics, during postyield loading in shock experiments. This behavior was termed “amorphization” in 2002 by Domnich et al.,¹⁶ who performed Raman spectral investigations on indented boron carbide and noticed the appearance of new peaks evolving near $1340cm−1$⁠, $1520cm−1$⁠, and $1810cm−1$⁠. With lack of additional data to confirm the validity, these new spectral features were labeled to emanate from graphite and carbon. These peaks, hence labeled D-band (⁠$1340cm−1$⁠) and G-band (⁠$1580cm−1$⁠), were assumed to appear due to the breakdown of boron carbide into carbon rings from broken C–B bonds. Contemporarily, in 2003, Chen et al.¹⁵ published TEM images of ballistically-impacted specimens of boron carbide. This was the first time that near-atomic-level imagery confirmed the presence of lattice disorder localized in tiny islands, similar to those in Fig. 1(a), stretching 1–2 nm in width and tens of nanometers in length, while, the material surrounding the disordered zones appeared overtly crystalline, with almost no change in the lattice structure. Since then, the investigation of the amorphization phenomenon in boron carbide continued with refinement in understanding the concept. In 2006, Yan et al.²⁵ studied the effect of temperature on Raman spectra of amorphized boron carbide, where it was found that annealing recovered the crystalline spectrum, decimating the so-called amorphous peaks. In 2007, Ghosh et al.¹⁹ proposed a rationale for amorphization based on visible and ultra-violet Raman spectra, photoluminescence, and FTIR spectroscopy data, suggesting the formation of carbon and boron clusters as part of the amorphization process, a proposition supported by DFT calculations of Gibbs free energy by Fanchini et al.¹⁰ Early investigations²⁶ discussed that amorphization might have some correlations with strain rate. However, it was later demonstrated that amorphization is, in fact, a rate-independent phenomenon,²⁰ and its intensity or extent depends on the applied pressure instead. In 2009, Yan et al.¹⁸ performed in situ Raman spectral measurements of compressed boron carbide using the diamond-anvil cell and concluded that boron carbide could withstand hydrostatic stresses up to 50 GPa, without any measurable structural breakdown, and it would rather amorphize (evident from Raman spectrum) when depressurized from high nonhydrostatic stresses and that the amorphization event occurred at $∼20GPa$ during the depressurization process. Subsequently, the structure of the amorphized zone beneath an indented region was determined by Subhash et al.²⁷ through successive material removal and polishing from the indented surface. Raman data collected from different depths were stacked to create a three-dimensional map of the amorphized zone. These studies were continued by Parsard and Subhash,²⁸ who concluded the Raman map to be independent of loading rate and that the amorphized state causes a local increase in volume, which leads to residual stresses in the surrounding matrix. These authors also found that in addition to the appearance of new peaks that primarily characterize amorphization, the dominant crystalline peak undergoes a shift as well [as shown in Fig. 1(b)], which has been related to residual stress. Postprocessing of TEM images through the generation of fast Fourier transform (FFT) diffractographs, published by Chen et al.,¹⁵ Reddy et al.,²⁹ and Zhao et al.,³⁰ has proved to be another evidence of amorphization. These studies show the presence of lattice order represented by regular dots, indicative of periodicity, while for amorphized regions, the analysis produces a diffused patch or halo, demonstrating lack of atomic-structural order. Recent advances on developing amorphization-resistant boron carbide⁹ highlight the effect of nanograining, which substantially reduces the propensity for amorphization in boron carbide.

Additional studies have been performed using the diamond-anvil test setup, notable among which are cases that did “not” claim or report amorphization. Fujii et al.³¹ used the diamond-anvil cell technique to study the effect of hydrostatic compression up to 126 GPa on structural changes through X-ray diffraction. While the X-ray diffraction profiles showed gradual peak-broadening and evolution of new peaks with pressure, authors claimed that amorphization was absent. No data on the unloading cycle were provided. Dera et al.³² studied X-ray diffraction patterns of single-crystal boron carbide for hydrostatic pressure up to 74 GPa, and while the main goal was to investigate the influence of compression on the nature of bonding of individual bond-types, no information was reported pertaining to unloading of sample, and no amorphization was reported to have occurred. Recently, works of Hushur et al.,³³ discussing optical properties, and compilation of experimental data from Werheit et al.,³⁴ focusing on the evolution of mode Grüneisen parameters, report that hydrostatic pressurization up to over 70 GPa leads to phase transitions without discernible amorphization, yet, in both these studies, depressurization behavior was not examined.

While the experimental investigations focused on the structure and spectral characteristics of amorphous zones, computational works focused on the mechanics of the amorphization process. In 2006, Fanchini and co-workers¹⁰ performed DFT-based Gibbs free-energy calculations to discover possible origins of amorphization in boron carbide. It was found that the polymorph (⁠$B12)CCC$ has the highest propensity among others to form a metastable structure that degenerates into alternating layers of graphite and all-boron icosahedra. It was also found that likelihood of achieving the segregated boron icosahedra and graphite phases increases with pressure and that other polymorphs, e.g., (⁠$B11Cp)CBC$ likely underwent distinct site-exchange steps among specific carbon and boron atoms, in order to achieve the configuration of segregated graphite phase. DFT studies of deformation by Aryal et al.²³ revealed bending of chains at high pressures to form new bonds between end-carbon atoms. This effect, where chain bending was identified to be of central importance, was found to cause an increase in the local volume and interpreted as the cause of amorphization. Few more DFT studies have shown that the application of compressive stress leads to chain bending,^14,35–37 demonstrating the chain to be susceptible to bending; potentially, the first step in the amorphization process. However, to date, no atomistic-level work has been able to validate a deterministic mechanism connecting chain bending to the alternating layered structure proposed by Fanchini et al.¹⁰ The chain bending proposition was further put to the test by DFT work of Betranhandy et al.,³⁸ who rather predicted that the end-carbon atoms of the C–B–C chain form a C–C bond, while the central boron atom gets released into the interstitial space. Raucoules et al.³⁹ performed an extended study on the stability of defects in the boron carbide system, concluding that the central boron atom in the C–B–C chain becomes susceptible to the detachment at pressures in the range 30–40 GPa, causing sudden volume drop. This process was attributed to be the driver for the loss in mechanical strength. There is another category of DFT-based stress-analysis of boron carbide which makes use of identifying slip systems with the weakest shear strength.⁴⁰ Central to these studies is the basic framework of constraining the crystal in a series of incremental shear deformations along different shear planes and evaluating the energy as a function of shear strain for each such plane.⁴¹ Originally developed for FCC aluminum and copper,⁴¹ this method has been used for several covalently bonded materials including diamond, cubic $BC2N$⁠,⁴² silicon, and germanium.⁴³ Even though new bonding is observed during the shearing process and higher shear strains led to complete breakage of the boron carbide structure, yet, these studies too have not been able to connect shear deformation to the alternating layered structure of Fanchini et al.¹⁰ DFT studies have also evaluated the Hugoniot curve for different polytypes³⁵ and compared them with experiments.⁴⁴ The comparison between experimental and computational predictions is very close for the most abundant polymorph, (⁠$B11Cp)CBC$⁠, in being continuous up to pressures of the order of 80 GPa. Other polytypes, including (⁠$B11Cp)CCB$ and (⁠$B12)CBC$⁠, which show discontinuities above 40 GPa in the pressure (⁠$P$⁠) vs volume ratio (⁠$V/V0$⁠) relationship, are speculated to participate in the amorphization process.³⁵ Moreover, discontinuity in the $P$ vs $V/V0$ dependence is interpreted as the onset of amorphization. The most relevant molecular dynamics (MD) work utilizes a recently developed ReaxFF potential for boron carbide,⁴⁵ where global shear deformations along the weakest slip system, namely, the $(0001)/⟨101¯0⟩$ plane were shown to produce twinning and the formation of a distinctive amorphous band. Though the ReaxFF potential used in Ref. 45 does not have a complete set of terms (e.g., absence of dihedral contributions from boron atoms), yet these simulations demonstrate shear-dominated deformation leading to microcracks within the amorphous band, which ultimately leads to failure.

In order to evaluate the pressure-dependence of Raman spectrum of boron carbide and decipher the origins of new peaks in the amorphized material, we embark on a systematic study involving the following steps:

1. Simulate the Raman spectrum of virgin boron carbide and compare with the experimental spectrum.

2. Simulate the evolution of spectral shifts due to applied hydrostatic pressure and compare them to those measured under diamond-anvil test results reported in the literature.

Though experimentally fabricated boron carbide is likely a mix of several closely related structures,¹¹ for simplicity, we present DFT-evaluated Raman spectra for the most dominant polymorph, $(B11Cp)CBC$⁠. To seek dependence of experimentally observed Raman peak shifts on deformation¹⁸ [see Fig. 2(a)] and relate them to applied pressure, we evaluated computational spectra for incremental levels of compressive hydrostatic pressure for $(B11Cp)CBC$⁠, which is shown in Fig. 2(b). Note that the DFT-evaluated spectrum for undeformed $(B11Cp)CBC$ [Fig. 2(b), labeled 0 GPa] depicts dominant peaks of experimental spectrum to be well captured. Jay et al.⁴⁶ were the first to evaluate Raman-active spectrum of pristine-crystalline boron carbide using DFT. Though most of their frequency set compare with our results [Fig. 2(b), 0 GPa case], we observe that Jay et al. did not calculate intensities, and their Raman spectrum did not possess three low wave-number peaks below $500cm−1$⁠. Our intensities were calculated for the temperature of 300 K and obtained using the method outlined in Refs. 47 and 48. For more details about the procedure of evaluating Raman spectra, the reader is directed to Sec. VI. DFT-evaluated Raman analysis was first made by Ivashchenko et al.⁴⁹ In that work, atomic models of boron carbide containing up to 135 atoms were annealed up to 4500 K temperature to produce amorphous material. DFT-evaluated Raman spectra of these systems, however, did not show much semblance to that of pressure-associated amorphized boron carbide and, essentially, none of the new spectral features of amorphized boron carbide (1340, 1520, and $1810cm−1$⁠) could be resolved. It is evident from Fig. 2(b) that high frequency modes (⁠$>500cm−1$⁠) undergo positive shifts, low-frequency modes (⁠$<500cm−1$⁠) undergo negative shifts, and the ones close to $500cm−1$ only show slight migration to the right with pressure. These DFT results have good agreement with several unique features observed in the experimental spectrum of hydrostatically compressed boron carbide in a diamond-anvil cell.^18,21 The specific features of interest, marked $f1$⁠, $f2$⁠, $f3$⁠, and $f4$ in Figs. 2(a) and 2(b), indicate the following features:

1. The low-intensity peak near $400cm−1$ (⁠$f1$⁠) evolves in intensity and moves to left with an increase in hydrostatic pressure.

2. The peaks near $480cm−1$ (⁠$f2$⁠) evolve gradually and move slightly to the right with increasing pressure.

3. At higher pressures, a new peak (⁠$f3$⁠) evolves at $1050cm−1$ and drifts to the right.

4. The peak at $1098cm−1$ (⁠$f4$⁠) increases in intensity and moves to the right.

The good comparisons between DFT-predictions of both overall Raman spectra [Fig. 2(a)] and pressure-dependence of Raman spectra with the experimental results for $(B11Cp)CBC$ [Fig. 2(b)] indicate a high degree of confidence in our approach. These results demonstrate the utility of quantum mechanical calculations in simulating the Raman spectra of boron carbide and quantitatively capturing the experimentally-observed trends. In the next steps, we take reference of DFT results to rationalize the unique behavior of experimentally obtained Raman spectrum of amorphized boron carbide, seen in Fig. 1(b). [Note that experimentally fabricated boron carbide is a mix of several polymorphs and contains varying degrees of defects. Hence, not all characteristics of $(B11Cp)CBC$ may exhaustively represent those of fabricated boron carbide.]

We recall the primary work of Fanchini et al.,¹⁰ and others²⁶ who embraced the idea of amorphization of boron carbide resulting in alternating layers of $α$-boron (⁠$α-B12$⁠) and graphite, and $(B12)CCC$ being the polymorph with highest susceptibility to amorphization of boron carbide. Adhering to that theory, in the neighborhood of the amorphization zone in boron carbide, we must expect a mixture of these three components: $α$-boron, $(B12)CCC,$ and graphite. The sequence of steps leading to amorphization can be summarized as⁵⁰ 

In order to verify if these reaction products indeed constitute the amorphous phase, and deepen our understanding of spectral evolution in boron carbide, we now simulate the Raman spectra of these components in stress-free states using the same approach [DFT and density functional perturbation theory (DFPT)] as we did for $(B11Cp)CBC$⁠. The simulated Raman spectra of these three materials and the experimental spectrum of amorphized boron carbide are shown in Fig. 3(a). We notice that simulated spectra do not conform to experimental spectral behavior of amorphized material. The spectra of $(B12)CCC$ and $α$-boron lie to the left of the prominent $1340cm−1$ feature. However, if $(B12)CCC$ and $α$-B$12$ are indeed products of amorphization (as predicted by Ref. 10), their spectra must rather exist toward the right to coincide with the amorphous band at $1340cm−1$⁠, hinting the components to likely exist in compressive states of stress. Also, the simulated spectrum of stress-free graphite does not occupy proximity to any of the spectral bands of amorphized boron carbide.

Motivated by the lack of overlap between the simulated spectral response of intermediate products of amorphization and the experimental spectrum of amorphized material, we incrementally subjected the unit cells of above three products to different levels of hydrostatic pressure, which caused their prominent Raman-active frequencies to right shift. In Fig. 3, we also show spectra at specific high hydrostatic pressure states for $(B12)CCC$⁠, $α−B12$⁠, and graphite at frequencies that closely match the experimental amorphous peaks. When all the spectral peaks are superposed (Fig. 3), interestingly, $(B12)CCC$ at 80 GPa hydrostatic pressure coincides with the most dominant peak near $1340cm−1$⁠, while another Raman-active peak of $(B12)CCC$⁠, situated around $1500cm−1$⁠, spans the spectral-neighborhood which has been referred to as the G-peak.¹⁶ Most interestingly, at the same pressure of 80 GPa, Raman peak of graphite is situated near $1800cm−1$⁠, the spectral “hump” which hitherto has never been resolved. Additionally, $α$-boron hydrostatically stressed at 105 GPa also occupies the $1340cm−1$ spectral zone. The pressure values quoted here are chosen based on the calculated peak matching the center of the broad amorphization peak. Given those conditions, is it correct to claim that the amorphous mix, composed of $αB12$⁠, $B12CCC,$ and graphite, exists in compressive pressure of the range 80–120 GPa?

The contemporary theory of the existence of intermediate products¹⁰ loses merit when we argue it from the perspective of our DFT results. If an amorphized island were to be under high compressive stresses (see Fig. 3, results for high-pressure case), where does the compressive stress originate from? If the surrounding crystalline matrix exerts such large stress on the islands that would result in the accumulation of equally high tensile stresses (up to 100 GPa) in the matrix, causing cracks in the immediate neighborhood of the amorphization region, leading to separation of the amorphous region from the surrounding crystalline matrix material. This will hold true for all other works which claim that amorphization results in the increase of density, e.g., Refs. 18 and 45, or reduction in the volume of the crystalline material that has amorphized. Additionally, subjecting a material mix containing graphite to high pressures may likely result in the formation of diamond, as recently shown by Blank et al.⁵¹ On the contrary, it must be noted that cracks or evidence of separation have never been found in the large volume of TEM images [e.g., Fig. 1(a)] of amorphized regions in the boron carbide literature nor have spectral characteristics of diamond ever reported in the amorphous spectrum of boron carbide. Thus, we prove that neither the amorphized region would be under compressive stresses nor possess higher density than the parent material.

Note that Raman spectra of $α−B12$⁠, $(B12)CCC,$ and graphite were calculated individually and those components are not mixed in any manner in the computational domain. The goal is to subject to test, whether the amorphous zone contains any constituents proposed by Fanchini et al.¹⁰ By the above analysis, we develop arguments that experimental Raman spectrum of boron carbide is not emanating from the three constituents: $α−B12$⁠, $(B12)CCC,$ and graphite.

Recalling that amorphization is induced in the presence of high pressure, we resort to high pressure compression of boron carbide using MD in order to observe the resulting temperature rise in the material. The strategy is to subject an MD model of boron carbide to different levels of volumetric compression, followed by complete relaxation, and observe the resulting structure. (Refer to Sec. VI for details of the MD setup.) Figures 4(b)–4(e) show structural features due to volumetric compression at $ΔVV0=15%$⁠, $20%$⁠, $25%,$ and $30%$ on the left and unloaded structure from that state on the right. [Note, the snapshots of structural features in Figs. 4(b)–4(e) are obtained from a 6 Å slice of the computational domain shown in Fig. 4(a).] We notice that as the imposed volumetric constraint is increased, temperature rises, manifested by high velocities of the atoms. At low volumetric compression (up to 20%, $T∼725K$⁠, and $P∼30GPa$⁠), the integrity of the structure is mostly maintained in both strained and released (unloaded) states, but at volumetric compression of about 25% (⁠$T∼1400K$ and $P∼40GPa$⁠), the structure loses most of its original order, having lost all its structural integrity even after unloading [see Fig. 4(d)]. When volumetric compression rises to about 30%, temperature and pressure reach 2100 K and 70 GPa, respectively. At this point, the atomic structure is quite different from that of the parent. When unloaded, the material sustains that disorder and stays amorphized. These studies conclude that high pressures of the order of 40–70 GPa can cause the onset of melting of the boron carbide structure and that melt, whether complete or partial, never regains its parent icosahedral configuration or its rhombohedral lattice, even after cooling. It is justified to refer to this high-temperature, high-pressure state of material as “melt,” because of two reasons: the loss of structural integrity and being at a temperature above the melting point. We argue that such conditions (⁠$P∼70GPa$⁠) occur in mechanical¹⁵ and laser shock experiments,³⁰ depressurization experiment,¹⁸ and indentation experiments,^9,16,27 where high pressures exist due to the application of mechanical loading.

Beyond the visual evidence of Fig. 4, radial pair distribution function (RPDF) is plotted for unloaded material states from different levels of volumetric strain, shown in Fig. 5(a). We observe that the RPDF of the crystalline state consists of unique peaks, of which the bonded ones belong to chain neighbors (at $r=1.42Å$⁠) and icosahedral neighbors (at $r=1.74Å$⁠). Populations pertaining to second neighbor distance appear clustered at $r=3$ Å. On application of incrementally higher levels of deformation, bonded populations reduce, as depicted by labels $R1$ and $R2$ in Fig. 5(a). We also observe the second neighbor populations to reduce, and in the process, a region between 2.2 Å and 2.7 Å (labeled by $R3$ in the inset) gets populated as shown by the arrow in Fig. 5. We notice an overall trend of reduction of peaks and leveling up of the RPDF with an increase in strain levels, confirming the incremental accumulation of amorphous state. Analyzing the drop in chain (⁠$R1$⁠) and icosahedral (⁠$R2$⁠) bonds, we conclude that in this simulation experiment, about one third of icosahedral bonds and about two thirds of chain bonds permanently break. We see a saturation of amorphization at peak strain levels of 35%, beyond which the RPDF remains unaltered. While such behavior might appear to limit amorphization in MD calculations, we argue that it is caused by the MD-methodology deployed in creating chemical bonds, being primarily dictated by spatial proximity. In reality, however, covalent bonding and rebonding events are likely more complex and depend on additional constraints.⁵² Therefore, the degree of amorphization predicted by MD is likely to be a lower bound.

Our results depicting RPDF of boron carbide possess similarity of bonded neighbors to those developed by Ivashchenko et al.⁴⁹ and Aryal et al.,²³ who also evaluated RPDF of amorphous models [Fig. 5(b)]. In both the above works, model systems were simulated using DFT and contained less than 200 atoms. Despite the wide difference in the number of atoms (18 000 vs 200), in essence, the amorphized state in each case consists of the main population level centered at $r=1.42$ Å and a second neighbor population at $r∼2.7Å$ [compare Figs. 5(a) and 5(b)]. A closer examination, however, reveals that RPDF for MD [Fig. 5(a)] demonstrates a greater degree of amorphization, as the second peak is much smaller than the first for MD as compared to the relative heights of the two peaks developed by DFT models [Fig. 5(b)]. Additionally, the third neighbor population for the case of MD is less prominent compared to DFT systems too. Ideally, complete amorphization would exhibit no peaks in RPDF.

Covalent bonding is strong and directional in nature; however, a system of unbonded atoms can be held together by weak nonbonded van der Waals type interactions, which can be attractive or repulsive, depending on separation between atoms. As part of the present theory, we derive closed-form analytical estimate of pressure associated with unbonded amorphized atoms, held together by van der Waals interactions, constrained to be packed in hexagonal-ordered arrangement and possessing the same density as that of the parent material (i.e., 2.52 g/cc for boron carbide). Conceptually, in this simplified model, the constraint of density limits the number of atoms per unit volume and the goal of placing them in regular arrangement predicts their spacing. Once the interatomic distance is determined geometrically, the degree of attraction or repulsion is assessed by the van der Waals potential. Interatomic force hence calculated from the interatomic separation can be further postprocessed using geometric considerations to yield the corresponding stress within such amorphous clusters. Details of such closed-form derivation of force and pressure of unbonded close packed atoms using nonbond van der Waals parameters are given in the supplementary material.

Such calculations yield interatomic separation of about 2 Å for the case of boron and carbon atoms packed at a density of 2.52 g/cc. That magnitude of interatomic separation gives rise to repulsive interactions, hence causing a pressure magnitude ranging between 150 (pure carbon) and 450 GPa (pure boron). In summary, if a given volume of boron carbide amorphizes without a change in volume, the amorphized atoms exerts high internal pressure on the boundary. This derivation assumes no motion of atoms, hence influences of temperature are not taken into account. In a real system, due to the thermal motion of atoms, the pressure would likely be higher than the above magnitudes.

There is currently a lack of understanding on what leads to bond breaking during amorphization, and what additional manifestations of amorphization are, especially to explain the presence of residual stresses in amorphized regions and how that influences the overall Raman spectrum of amorphized material. Results in Secs. III A–III D necessitate a combined explanation of amorphization in boron carbide that can relate experimental evidence of TEM, namely, the presence of amorphous material in islands to the observation of new features in the experimental Raman spectrum. Our novel approach consists of elucidating the pressure-induced amorphization process to be thermodynamic phenomena, beyond mere breaking of bonds. Using arguments of temperature rise when subjected to high pressure, we rationalize that local changes in crystalline morphology are likely to occur in confined pockets of material near inhomogeneities and defects, manifesting in local melting events. The melt, when cooled rapidly, is less likely to bond as parent crystallite and contains unbonded atoms, which exert strong repulsive forces on each other and the surrounding matrix, thereby causing residual stresses in the crystalline matrix. Finally, we demonstrate that Raman spectra collected from the crystalline matrix under residual stresses provide a good comparison with experimental Raman spectrum of amorphized material.

An inherent weakness in reported contemporary computational models is that they all refer to the formation of intermediate crystalline phases: whether bent-chain morphologies^23,36 or having proposed a mix of $(B12)CCC$⁠, $α−B12$⁠, and graphite,¹⁰ lack evidence of their widespread dispersion in experimental TEM samples. Such crystalline by-products of the so called amorphization process must rather be termed “phases” and their appearance must instead be termed “phase transformation,” not amorphization. While phase transformation refers to one deterministic structure transforming to another deterministic structure, amorphization, as observed in TEM, refers to the destruction of the deterministic structure, leading to a lack of long-range order. Some recent atomistic studies link bond-breaking events to be connected to amorphization.^40,45 Application of high shear loads along directions of weakest shear strength will likely disrupt parent bonding, which is the main strategy of inducing amorphization in these studies. Mechanisms of bond breaking in amorphization (as we discuss in Secs. IV A–IV C) is different from damage through mechanical means and must hence be treated differently. Since the early shock studies by Grady,²⁴ which were continued further by Vogler et al.,⁵³ there has been a tendency in the community to seek “phase transformation” behavior of boron carbide at high pressure, something which was never conclusively determined.

Past works that have shown neither amorphization nor presence of any drastic instabilities for pressures up to and over 70 GPa^31–34 have all deployed quasistatic hydrostatic deformation. It is important to note that in these studies, the dominance of shear, which is a major requirement for amorpization, is not well incorporated, because of the hydrostatic conditions at macroscopic level. In addition, from the perspective of atomic-level understanding, we can further notice that quasistatic loading ensures equilibration of pressure and temperature across length and time scales, hence even if any amorphization nucleated, it would very soon dissipate, bringing back the primary crystal structure.

In view of lack of phase transformation, and clear evidence of amorphization, new directions toward understanding the mechanical behavior of boron carbide at high pressures are sought for. To achieve that goal, the motivation of the present work is to seek a new approach to rationalize deformation behavior of boron carbide that takes material from deterministic structure to undeteministic configuration and be consistent with postamorphization conditions, such as explanation of residual pressure in surrounding crystalline matrix, as well as appearance of new features in the Raman spectrum.

Now that we have sufficient understanding of Raman spectra of boron carbide under pressure, along with a link between pressure and solid and melt phases, let us augment our understanding with thermodynamics of the amorphization process itself. It is known that high pressure not only raises the homogeneous temperature [see $Thomo$ in Fig. 6(a)], but also causes reduction of melting point [see sloped line labeled $Tm$ in Fig. 6(a)]³⁰ and at those pressures, regions close to inhomogeneieties convert a significant portion of work into heat, thereby raising the temperature of those pockets close to melting. Hence, in a material subjected to macroscopic applied stresses of the order of 50–60 GPa, defects or other inhomogeneieties present in the material act as sources for stress-concentrations and may cause the pressures to rise much higher, causing local melting [Fig. 6(a)]. Given the loading timescales operating in these experiments (namely, shock, indentation, and depressurization in the diamond anvil), the melt, which is formed at high pressures, possesses high atomic velocities. We conjecture, that the melt does not get enough time to cool back to the parent configuration and hence ends up as amorphous solid. As shown in Fig. 6(a), the regions where amorphization occurs, experience stress well above the average HEL conditions. This is because, unlike plasticity, amorphization does not occur homogeneously throughout the deformed zones, but in sporadic islands [see Fig. 1(a)] near defects which experience high stress concentration compared to homogeneous macroscopic stress. This stress enhancement (concentration) can cause a significant rise in temperature as the applied macroscopic stress increases. Hence, the presence of defects reduces the threshold stress required to cause local melting. This effect is represented by arrows in Fig. 6(a), highlighting that in the presence of defects that possess stress intensity factors greater than one, HEL is reached locally much earlier. By the same effect, local melting events occur at smaller global stress. Due to this reasoning, it is easier to understand that indentation experiments demonstrate amorphization, because even though global forces are low, they are sufficiently significant to cause melting and hence amorphization.

Note that parameter $β$ represents a fraction of plastic work done by shear converted into heat, and that $β$ increases with the magnitude of shear. Quasistatic hydrostatic loading is isothermal in nature and likely represents negligible $β$⁠. In our MD simulations of hydrostatic loading, time scale of loading is very small (order of picoseconds) and represents an adiabatic condition. In a given experiment of subjecting a sample of boron carbide under high pressures, where both amorphous and crystalline phases exist, temperature rise pertaining to the entire macroscopic sample is not as drastic as that in the localized amorphous zone. Our MD results pertain to material behavior represented by $Tband$ of Zhao et al.³⁰ [refer Fig. 6(a)] and correspond to conditions within the amorphous zone. Comparison with temperature rise presented in Fig. 6(a), we infer that our MD simulations pertain to $β$ value of about 0.6. Extending our understanding of the interplay between shear and temperature causing irreversible structural damage in boron carbide, we emphasize that diamond anvil tests of hydrostatic loading with quasistatic loading rate, especially the work of Fujii et al.,³¹ Dera et al.,³² Hushur et al.,³³ and Werheit et al.,³⁴ represent very low $β$ and are unlikely to produce permanent amorphization, because temperature rise associated with such deformation is relatively low and merely governed by $Thomo$⁠.

Let us also look at effects of rates of loading. Consider a block of boron carbide compressed quasistatically and hydrostatically from ambient conditions to some high-pressure value above the HEL (say 50 GPa). Thereafter, consider the following three separate situations:

• • Case A: external loading is maintained at the same pressure (of 50 GPa),

• • Case B: external loading is continued by increasing pressure, and

• • Case C: external loading is rapidly withdrawn.

Before being subjected to any of the above cases, the compressed block will likely contain several tiny islands dispersed within the volume, whose temperature is above the melting point. Transfer of heat from molten islands to surrounding solid material establishes a phase equilibrium between solid and molten counterparts. Atoms in molten islands have much higher kinetic energy than those in solid phase. Such atoms, especially those with broken bonds, can eventually rebond with the correct neighbors (recrystallization) in Case A, by consequence of equilibrium, given sufficient time. On the other hand, Case B consists of interplay of two drivers: temperature rise (leading to melt) and recrystallization (leading to resolidifcation), which depends strongly on the rate of loading. If the loading rate is sufficiently high, or there is a predominance of shear, pronounced melting will lead to structural collapse, while for lower loading rates, equilibrium and recrystallization will likely support structural integrity. This explains why it is possible to load boron carbide quasistatically to over 120 GPa without structural collapse.^31–34 In Case C, it is quite likely that the atoms in molten state do not get sufficient opportunity to recrystallize, because rapid unloading leads to rapid temperature drop, causing enhanced solidifcation, due to which, atoms have reduced mobility to reach out to potential neighbors to recrystallize. Thus, Case C likely produces amorphous zones, constituting unbonded atoms in solid state.

At this point, it is also important to recognize the distinction of processes in molten zones of boron carbide at high temperatures compared to those of high-temperature molten phases of metals. In metallic systems, during cooling from high temperature, rebonding is usually straighforward, because atoms in molten state seek satisfaction of parental coordination. Once sufficient neighbors aggregate in the necessary proximity for rebonding, coordination-driven processes give rise to recrystallization. However, in covalent materials such as boron carbide, during solidification, atoms in molten state seek not only parent neighbors, but also their parental “orientations,” an added constraint that makes recrystallization not as efficient as in metals. Hence, the tendency of boron carbide to exist in amorphized states post solidification. Because of broken bonds, these amorphized zones have extremely low structural integrity even after solidification and are extremely susceptible to failure.

Referring to our analytical derivation of pressure due to unbonded atoms, it is now clear that a small domain of amorphized material consisting of unbonded atoms with repulsive forces, exerts high pressure on the surrounding crystalline medium. If the interface between the amorphized and crystalline regions are flexible, the confined amorphous system will likely expand (dilate) to alleviate the high pressure, thereby exerting compressive stresses in the surrounding medium. The equilibrium state of the resulting system is an expanded island surrounded by compressed matrix material, likely having continuity of tractions at the interface.

From the reference frame of the surrounding crystalline matrix, this effect is similar to volume expansion of the amorphous island, leading to a reduction in its density. Previous experimental works (Parsard et al.²²) have argued the presence of residual stresses in the matrix material due to volume expansion of amorphous islands. However, those studies never reasoned the source for volume increase. The current theory rationalizes these results based on thermodynamics and nanoscale force considerations.

The present theory evolves from the basic assumption that amorphized regions are constituted of unbonded atoms which, having severed from the parent structure, possess no chemical bonds to produce any spectral signature. We argue that the so-called “amorphized spectrum” emanates from crystalline regions surrounding the amorphized zones. In this section, we discuss the impact of residual compressive stresses within the matrix material to the overall Raman spectrum observed in an experiment. Currently, the laser beam used in Raman instruments is limited in diameter to $1μm,$ and, hence, cannot effectively characterize and resolve nanoscale amorphous material regimes, which are nanometer size, as shown in Fig. 6(b). The same figure also shows the proposed neighborhood of amorphized zones, consisting of amorphous atoms, interface, and vast crystalline material surrounding the amorphous atoms.

Furthermore, the amorphous material exerts pressure on interface, compressing it from inside, causing residual pressure in the matrix material. Hence, the presence of nanosized amorphous material leads to superposed Raman spectra emanating from a spatial distribution of material pressurized from inside. That pressurized segment (composed of superposition of spectra pertaining to different stress states) has its own spectral contribution, while far away from that zone, the less-affected crystalline segment contributes its own pristine-crystalline spectrum. Thus, the overall spectrum contains a mix of (a) spectrum of pristine, unstressed material and (b) mix of spectra of crystalline material stressed at different stress levels, giving rise to features at 1340, 1520, and $1810cm−1$⁠. This approach is a new perspective of analyzing Raman spectra of material systems with submicrometer crystallographic entities, which cannot be resolved due to limitations of the present-day spectral resolution of the instruments. Additional aspects of this model are described below.

The compressive stress developed in the matrix material, which might be of high magnitude near the boundary between amorphous island and the matrix, decays with distance away from the amorphous zone. Additionally, the stress state in the matrix material is likely to be hydrostatic. Hence, we argue that the overall Raman spectrum of the amorphized material is likely to be a superposition of Raman spectra of matrix material stressed to different levels. Figure 6(c) shows superposition of Raman spectra of four different stress levels: 0, 5, 60, and 90 GPa. (Note, a low-frequency mode situated near $200cm−1$ is not shown because it becomes unstable at high pressure. This mode does not recover upon unloading and continues to stay unstable, and hence does not contribute to the overall spectrum.) Incidentally, many salient features of the experimental Raman spectrum of amorphized material seem to be captured, including the dominant $1340cm−1$ peak, the so-called graphitic peak at $1520cm−1$⁠, as well as the “unexplained” hump near $1800cm−1$⁠. We further argue that the intensities of these features would vary in proportion to the amount of volume of matrix material stressed. With an increase in the number density of amorphous islands, more matrix volume would get stressed, leading to higher intensities of the new spectral features. Thus, according to this theory, all features of the Raman spectra of amorphized boron carbide originate from differently stressed regions of the matrix material which is crystalline. It should be reiterated that the above computational spectra were obtained by considering only the dominant polymorph of boron carbide, namely, $(B11Cp)CBC$⁠. However, any fabricated boron carbide has potentially numerous polymorphs¹³ at lower volume fractions and hence some discrepancy is expected.

According to the current theory, upon application of high pressure, temperature rises up to close to melting point in localized regions which are surrounded by crystalline matrix material. Upon removal of loads (and impending solidification of melt), these zones undergo mismatch at their interfaces and the resulting amorphous material exerts compressive stresses on the surrounding matrix. Other efforts have merely hypothesized the source of these residual stresses to be connected with likelihood of local volume increase.²⁸ The present theory provides validation to that claim, as fundamentally, the process of bond breaking at high temperatures and solidification thereafter introduces lattice disorder, which is higher in energy than original undeformed crystal structure. Constrained by ordered crystalline material from outside, the amorphization zone exerts pressure from within, which is the reason for its locked state of stress. It is also important to note that even though the amorphization neighborhood possesses high residual stresses of the order of 80–120 GPa and HEL being close to 22 GPa [as per Fig. 6(a)^,30], the material would not continue to amorphize, because conditions of even higher stresses [see circled region of Fig. 6(a)] are required to initiate further amorphization.

One consequence of locked state of stress is that if stresses in the range 60–90 GPa exist in crystalline portion, the crystalline material will likely not re-amorphize. Based on the rationale of amorphization discussed in this paper, the amorphized region expands, exerting compressive pressure on surrounding crystalline material and that steady-state stress state in crystalline material is hydrostatic in nature. Studies made on boron carbide pertaining to quasistatic hydrostatic loading on diamond-anvil setup have shown^31–34 (as previously discussed in this paper) that this material can tolerate high hydrostatic pressure conditions up to over 120 GPa. Thus, the locked state of stress is stable in nature, and high residual stresses will not lead to further amorphization.

It is important to refer to the work of Yan et al.²⁵ who demonstrated experimentally, that the transition from crystalline to amorphous phases in pressure-induced amorphization can be reversed by heating the sample to over 800 K followed by cooling. This feature can be well explained by the present model, as we can readily judge, that the amorphized atoms, which separated from the parent structure and could not get back to their original locations, upon heating, are likely to absorb additional energy to increase their probability of restoring either their original neighbors or similar neighbors. Same response would be achieved if the cooling cycle after amorphization were to happen slowly. The present model also explains the small right shift of the dominant crystalline peak of $1080cm−1$ in the experimental Raman spectrum of amorphized specimen. This likely arises due to relatively lower levels of compressive stresses (about 1–5 GPa) in matrix material far from the amorphized zone.

The present work highlights that pressure-induced amorphization whether in indentation, shock, or diamond-anvil experiments consists of localized nanoscale melting events. Upon quick cooling, due to unloading of externally applied loads, the melt does not get sufficient opportunity to rebond like the parent structure, giving rise to residual stresses through unbonded atoms, which in turn alters the overall Raman spectrum. The drawback of the present model is that we can neither predict the degree of amorphization, nor the size of the amorphous zone. Perhaps, a large-scale MD simulation consisting of an embedded amorphous zone inside a crystalline matrix might be a good direction. That is, the most intuitive procedure to generate the final spectrum. However, there are a few limitations pertaining to the availability of software platforms. While Raman spectra are evaluated accurately from DFT, it is largely limited by size of the computational domain. Similarly, even though force field based methods can accommodate a larger number of atoms, yet their predictions on Raman spectra have not been established yet. We also need better force-field representation of materials like boron carbide, which contain accurate switching criteria between bonding and nonbonding interactions to model processes like amorphization more realistically at all length scales. Our model of deciphering the overall Raman spectrum of amorphized boron carbide assumes all spectral contributions to originate from matrix (crystalline) material, which is different from the current explanation of new spectral features to emanate solely from amorphized zones.^16,20,26 Lastly, distribution of residual stress is unknown and, therefore, the correct stress state of the crystalline medium is not known. Perhaps an Eshelby-like development⁵⁴ for this specific problem could also be attempted so that the overall Raman spectra could be accurately predicted.

Vacancies, like other inhomogeneities are likely to play a key role in the amorphization process. However, theoretical studies on their role in amorphization are not straightforward, and very few studies have been undertaken to describe their role in amorphization.^38,39 Vacancies give rise to stress concentration by causing local shrinkage,³⁹ thereby reducing the amorphization-threshold of macroscopic compressive pressure [depicted by arrows in Fig. 6(a)]. Moreover, investigating stability of boron carbide structures having vacancies is a challenging problem, making evaluation of theoretical Raman spectra and initial minimization of structures to do MD-level analysis harder. Though treatment of such topics is outside the scope of the present work, they are excellent directions for future.

In this paper, a new rationale has been presented to better understand the origin and manifestation of pressure-induced amorphization in boron carbide, a commonly studied icosahedral boron-rich ceramic. Pressure-induced amorphization in boron carbide has been observed in indentation, mechanical shock, and diamond-anvil experiments and possesses common TEM imagery and Raman spectra. Since its discovery, over two decades ago, this phenomenon has been shrouded in mystery, especially in (a) its occurrence as dispersed nanosized highly localized islands and (b) the appearance of new features in Raman spectra (at 1340, 1520, and $1810cm−1$⁠). Hitherto, it was believed that these new spectral features might emanate from free carbon, boron, and $(B12)CCC$ present in the amorphization zone. Using Raman spectra derived from DFT, we demonstrated that these components must be under high compressive pressure to produce spectral resemblance consistent with experimental observations. We also proved that such theories have an intrinsic flaw that amorphization zone under compressive pressure would lead to tensile stresses in the surrounding matrix, thereby producing cracks around the amorphization zones, which have never been observed experimentally and even formed diamond, which has never been observed in Raman spectra of amorphized boron carbide.

In developing our new theory, we have simulated Raman spectra using DFT for different stress levels and performed MD simulations of volumetric compression using the reactive force-field ReaxFF. The new rationale follows a thermodynamic basis, in that, amorphization occurs near zones of high stress concentration, reaching melting temperature due to applied pressure. Our MD simulations confirm that boron carbide melts at high compressive pressures of about 70 GPa. These simulations also show that when the applied pressure is released, the melt does not recrystallize back to the original structure and consequently remains amorphous. Using nonbonded interactions from the force field, we demonstrate that such amorphized atoms exert repulsive forces on each other and when present in a nanoscale confined region, exert compressive stresses of the order of 150 GPa on the surrounding crystalline matrix, which manifest as residual pressure in the matrix material. Crystalline matrix subjected to high residual pressure impacts Raman spectrum of crystalline region that manifests as new spectral signature unique to amorphized state. DFT-simulated Raman spectra for pressures between 60 and 80 GPa possessed same features as that of the experimental Raman spectra of amorphous boron carbide, thus proposing a new explanation for the origin of the mysterious new spectral features at 1340, 1520, and $1810cm−1$⁠. The new rationale of the appearance of new Raman peaks is summarized in the following:

Pressure-induced amorphization in boron carbide is caused in small nanosized pockets in regions near inhomogenities and defects under conditions of high pressure concomitant with shear. Amorphization occurs owing to rise in temperature above the melting point, caused due to high pressure, and is hence a thermodynamic process. When cooled quickly, the amorphized material fails to rebond covalently as parent crystalline material, and its nonbonded atoms exert high pressure on the surrounding crystalline matrix due to consequent volume mismatch and confinement. The crystalline material thus gets subjected to compressive residual stresses near amorphous zones and the state of stress gets locked due to the volume mismatch in the amorphous zone. Raman spectrum is mainly the overlap of the spectrum of differently stressed crystalline regions as a consequence of this locked stress state.

Our approach demonstrates the importance of multiscale considerations in evaluating material response, especially in circumstances where limits in experimental resolution cannot distinctly unravel the structure, composition and response of the material. This methodology is also important when spectral signature of an amorphized compound lies in close neighborhood of spectral signatures of its elemental components. In pursuit of ultrahardness, this study also outlines the necessity of discovering and rationalizing material behavior at extremely high pressures. It is important to identify the link between mechanisms originating at the atomic scale and their gross manifestation at the meso- and macroscales to develop a true understanding of properties of ultrahard materials.

Hot-pressed boron carbide was obtained from Ceradyne Inc., USA. Indentation was performed using the Tukon® 2100B Wilson Instruments hardness tester. Raman spectral experiments were carried out using the Renishaw® InVia Raman microscope along with the 532 nm (Si) laser.

Simulated Raman spectra for boron carbide was obtained using DFT and DFPT computations on fully periodic rhombohedral unit cell of boron carbide (15-atom). Spectra for pristine and hydrostatically compressed material were hence obtained. We utilized the open source software ABINIT,⁵⁵ using plane wave basis sets and norm-conserving pseudopotentials with the Troullier-Martins scheme and Pulay mixing. The exchange-correlation functional is evaluated by the Teter-Pade local-density approximation (LDA), and a $4×4×4$ Monkhorst k-point mesh and a cutoff energy of 680 eV are utilized for all simulations. All atomistic systems were first subjected to a conjugate-gradient energy minimization which is accurate up to $10−13eV$⁠. After the energy minimization, we recorded crystal lattice-parameters, which were in good agreement with previous DFT studies,^11,40 as well as experiments.⁵⁶ The energy-minimized structure was then subjected to the DFPT module of ABINIT⁵⁷ from which crystal vibration characteristics and Raman spectra were evaluated through postprocessing vibrational information as detailed in Refs. 47 and 48, using laser frequency of 532 nm and temperature 300 K. This methodology of obtaining Raman spectra of crystals is consistent with efforts of a larger community called WURM.⁵⁸ In presenting these results (both DFT and experimental), the spectra are normalized by the intensity of the most dominant peak. Spectral information output from ABINIT consists of natural frequencies and intensities but not peak width.

The MD simulation cell of boron carbide is a rectangular domain, prepared from the original 15-atom unit rhombohedral cell of the $(B11Cp)CBC$ polymorph. There are 18 000 atoms in the MD block, and it is periodic in all directions. Volumetric compression was carried out by reducing the size of this computational domain through incremental volumetric strain. This was implemented using the “fix deform” module of LAMMPS through which all the three box dimensions were incrementally shrunk by equal amount of strain in each step. Pressure and temperature were monitored during the compression process. At each level of volumetric strain of 15%, 20%, 25%, and 30%, the system was unloaded to restore it back to temperature and pressure conditions of 300 K and 1 atm, respectively, using the constant temperature constant pressure ensemble, also known as the NPT ensemble. The intent was to observe the residual structure after subjecting it to the above volumetric strain and allowing it to relax to original temperature and pressure conditions, similar to a typical indentation experiment where pressure is applied and removed. Note that the relaxation process was performed differently than the loading step (NPT run equilibrating to $T=300K$ and $P=1atm$⁠).

Internal pressure associated with unbonded atoms, confined to exist with the same density as boron carbide, i.e., 2.52 g/cc, is derived using nonbond force field parameters and presented as the supplementary material.

The support provided by Army Research Office under Grant Nos. ARO-W911NF-14-1-0230 and ARO-W911NF-18-1-0040 were highly acknowledged. The authors thank support from the Department of Defense pertaining to the acquisition of an EMCCD Raman spectrometer via Grant No. ARO-W911NF-16-1-0180. DFT- and MD-simulations discussed in the present work have used computational resources from the Extreme Science and Engineering Discovery Environment (XSEDE) program, which is supported by the National Science Foundation Grant No. ACI-1548562. Resources from research allocation No. TG-MSS160016 were used. The authors also thank Matthew DeVries for preparing the rectangular atomistic domain for MD simulation of boron carbide.