Profound softening and shear-induced melting of diamond under extreme conditions: An ab-initio molecular dynamics study

Carbon 155 (2019) 361e368

Contents lists available at ScienceDirect

Carbon journal homepage: www.elsevier.com/locate/carbon

Profound softening and shear-induced melting of diamond at extreme conditions: An ab-initio molecular dynamics study Libin Wen a, b, Hao Wu a, c, Hong Sun a, b, *, Changfeng Chen d, ** a School of Physics and Astronomy, Key Laboratory of Artificial Structures and Quantum Control (Ministry of Education), Shanghai Jiao Tong University, Shanghai, 200240, China b Collaborative Innovation Center of Advanced Microstructures, Nanjing, 210093, China c School of Physics and Electronics, Hunan Normal University, Changsha, Hunan, 410081, China d Department of Physics and Astronomy, University of Nevada, Las Vegas, NV, 89154, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 10 June 2019 Received in revised form 28 July 2019 Accepted 26 August 2019 Available online 28 August 2019

Diamond retains original crystal form and large elastic stiffness up to very high temperatures, and strengthens considerably under high pressures. These remarkable characters dictate its behaviors under extreme loadings such as those in laser-heated diamond anvil cells (LHDACs) or giant planetary interiors. Here we report surprising findings from ab initio molecular dynamics simulations revealing that diamond's mechanical strengths reduce precipitously at elevated temperatures, by 50% at 3000 K, despite that its elastic parameters remain little changed. Our results also unravel an extraordinary shear-induced reduction of melting temperature of diamond under ð111Þ½112 shear strains by as much as 1150 K, greatly altering its pressure-temperature (pT) phase diagram when shear deformations exist. These new benchmarks offer crucial insights for elucidating extreme mechanics of diamond at high p-T conditions, advancing fundamental knowledge with major implications for LHDAC operation and planetary modeling. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction The structural and mechanical properties of diamond at highpressure and high-temperature (HPHT) conditions are of great interest both for fundamental significance in understanding ultimate material characters and for practical importance in applications involving diamond as a major component. A conspicuous case is the use of diamond as the key component in laser-heated diamond anvil cell (LHDAC) [1]. One important area of application of LHDAC is new materials synthesis involving HPHT in the range of 10  100 GPa and 1000  3000 K [25]. Another significant area of application is geophysics, where LHDAC is the only instrument that can maintain steady HPHT in the range of 100  300 GPa and 2000  7000 K to study the structure as well as electrical and thermal conductivities of minerals likely existing in Earth's mantle or core

* Corresponding author. School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China. ** Corresponding author. Department of Physics and Astronomy, University of Nevada, Las Vegas, Nevada 89154, USA. E-mail addresses: [email protected] (H. Sun), [email protected] (C. Chen). https://doi.org/10.1016/j.carbon.2019.08.079 0008-6223/© 2019 Elsevier Ltd. All rights reserved.

regions [6e9]. Moreover, the structural stabilities and mechanical strengths of diamond also become crucial in astrophysics, where diamond has been found in meteorites reaching earth [10e12] and postulated to exist at extreme HPHT in the range of 300  600 GPa and 3000  7000 K, inside carbon-rich giant planets, such as Neptune and Uranus [13], which often requires experimental tests using shock wave compression [14e16]. Experimentally, the elastic and plastic mechanical properties of diamond at high temperatures have been extensively studied, and the results show that the elastic constants of diamond change slowly with temperature up to about 1600 K [17], while its plastic mechanical properties, such as hardness, decrease quickly as the temperature increases above 1000 K [18-20]. At extremely high pressures ( > 300 GPa), most studies mainly focus on the equation of state and pressure-temperature (p T) phase diagram of diamond [21e23]. It has been predicted that diamond remains stable in the core regions of Neptune and Uranus where pressure and temperature reach 600 GPa and 7000 K [13e16]. However, the mechanism remains to be elucidated on how mechanical strengths of diamond change in such an extreme HPHT environment, especially how shear stresses would affect diamond's phase stability and melting temperature. It is noted that shear-induced melting

362

L. Wen et al. / Carbon 155 (2019) 361e368

has been reported for a variety of materials, such as colloidal suspension and dusty plasma systems, and adopted as an explanation to the anomalies in melting behaviors of some metals, such as Ta etc. [24e27]. It has long been a great challenge to establish an atomistic understanding about the mechanical strengths, phase stability and melting temperature of diamond under shear stresses that commonly exist in many realistic HPHT environments, such as those in LHDACs and inside meteorites and carbon-rich giant planets. Currently, theoretical studies on stress-induced melting are mostly based on classical molecular dynamics (MD) simulations, which employ classical interatomic potentials to describe atomic motions, and the reliability of the results heavily depends on the quality of the potential functions. Previous studies using classical MD show important and convincing findings on single-element crystals like Cu [28], Al [28e29], and systems with simple interatomic potentials (such as the Lennard-Jones potential) [30]. However, it remains highly challenging to develop accurate and universal potential functions for situations where breaking and recombination of atomic bonds frequently occur, especially in complex compounds, or large density variations exist [22] in the simulation process. Such shortcomings can be overcome by using ab initio molecular dynamics (AIMD) approaches, where forces acting on atoms are determined by more accurate quantum mechanical calculations. In this work, we present findings from AIMD calculations of mechanical strengths and deformation microstructures of diamond at high temperatures, focusing on its phase stability and melting temperature under shear stresses at extreme HPHT. Our results show that under ambient pressure, mechanical strengths of diamond are reduced by about 50% as the temperature reaches 3000 K, while its elastic constants remain almost unchanged, which is consistent with the experimental results [17e20]. At high pressures (100  600 GPa), we adopt the small-size coexistence method developed by Hong et al. [31,32] to determine the melting temperature Tm of diamond under hydrostatic pressures, followed by an examination on phase stability of large-shear strained supercell at HPHT. The obtained results show that under large shear strains, carbon liquid state becomes more stable than shear strained diamond crystal state at a shear-induced melting temperature that is much lower than the melting temperature of the unstrained diamond Tm , with the temperature reduction as large as 1150 K. These findings establish new benchmarks and provide crucial insights for understanding extreme mechanics of diamond under extreme HPHT conditions. 2. Methods of Calculations Our calculations were performed using the AIMD code implemented in the Vienna Ab initio Simulation Package (VASP) [33e35], adopting the projector augmented wave potentials [36] and the generalized gradient approximation with the exchange-correlation functional by Perdew, Burke and Ernzerhof [37] for the calculations of the electronic states and forces on atoms, where a plane-wave basis set with a 500 eV energy cutoff and 2 2  2 MonkhorstPack [38] k-point grid were used for a cubic diamond supercell of 64 carbon atoms with periodic boundary conditions (see Fig. 1(a)). To obtain the equilibrium structure at a given temperature and pressure ðT;PÞ, we start from the equilibrium structure at zero (T;P), and then increase (T; P) of the previously obtained equilibrium structures step by step, with T ¼ 300 K, 1000 K, 2000 K, … and P ¼ 0 GPa, 100 GPa, 200 GPa, 300 GPa, …, to the target (T;P). The diamond equilibrium structure at each (T, P) is determined by the Parrinello-Rahman (NpT) AIMD [39] simulation runs for at least 3000 steps (6 ps) with the temperature controlled by the Langevin

thermostat [40], until the equilibrium of the simulation is reached, for which we can monitor the averaged stress components and structure parameters of the supercell. Although the temperature T increases directly from the previous equilibrium structures step by step, the effect of the temperature gradient to the outcome structure is minimized by the full relaxation at each equilibrium structure, where both the convergence and equilibrium are verified. It should be noted that, however, for system containing complex carbon phases or structures under fast heating or quenching processes, the temperature ramping rate may extensively affect the final outcomes [41], and will need more careful considerations. For the calculation of the stress-strain relation and loading path under compression and shear stains, we follow the procedures of our previously developed quasi-static DFT approach at zero temperature [42e44]. Starting from the average equilibrium structure at a given pressure and temperature, we continuously deform the supercell to achieve the target shear strain (εzx ) or compression strain (εzz ) at a constant strain rate ε_ ðtÞ ¼ 0:001ps1 (see discussion below). The motion of the atoms and the shape of the supercell are determined by NpT AIMD [39] with a time step of 2 fs. During the straining process, the applied strain component, taking shear strain εzx as an example, is controlled by the strain rate ε_ ðtÞ and AIMD time steps, while the other five independent strain components of and all the atoms in the supercell are free to move according to the NpT AIMD simulation process, which satisfies the following conditions: (i) The average (say over 2000 AIMD steps) of the uniaxial stress components sxx , syy and szz equals the hydrostatic pressure p; (ii) The average of sxy or syz is nearly zero (see Fig. 1(e) where p ¼ 0 GPa). The temperature during the AIMD stress-strain simulation is kept as that of the equilibrium structure by the Langevin thermostat [40]. The (average) shear stress szx , the shape of the supercell, and the motions of all the atoms in the shear deformed diamond structure at a given hydrostatic pressure and temperature as a function of the applied shear strain εzx , which is linearly proportional to the AIMD time steps, are fully determined by the AIMD simulation. To calculate the melting temperature or check the thermodynamic stability of the (strained) diamond structure by AIMD, we use the small-size coexistence method developed by Hong et al. [31,32], which has been applied successfully to predict melting temperatures of various materials [45,46]. Starting from the AIMD averaged (over 2000 steps) diamond structure under given pressure and temperature with or without applied (shear or compression) strains, the supercell is doubled along the ½001 direction, then the solid-liquid coexistence configurations are prepared by heating and melting half of the doubled supercell, while keeping the other half of the atoms fixed in their initial solid positions. Using these coexistence configurations as initial structures, at least 3000 steps of the NpT AIMD simulations described above are carried out at given temperature to monitor the evolution until pure final phases are reached, either solid or liquid. Melting temperature can then be deduced from the ratio of the number of final liquid states (NLiquid ) to that of solid states (NSolid ) at various temperatures, a≡NLiquid =ðNLiquid þ NSolid Þ. If a > 0:5, the liquid state is more stable than the initial (strained) solid state. The melting temperature is considered reached when a ¼ 0:5 [31,32]. 3. Results and discussions In Fig. 1(b), we plot calculated elastic constants, bulk and shear moduli of diamond at different temperatures under zero pressure in comparison with experimentally measured elastic constants. Our results agree well with those of experiments [17] and previous calculations [47,48], which show less than 15% changes in the elastic properties of diamond at temperatures up to 3000 K. In

L. Wen et al. / Carbon 155 (2019) 361e368

363

Fig. 1. Mechanical properties of diamond at different temperatures under 0 GPa. (a) A schematic illustration of the shear and compression directions of the diamond supercell used in the calculations. (b) The elastic constants, bulk moduli and shear moduli. Data at 0 K is obtained from DFT calculations, while data at finite temperatures are obtained using AIMD. Experimentally measured data are also plotted for comparison [17]. (c,d) The average stress-strain relations of diamond at different temperatures under shear deformation along the ð111Þ½112 slip direction and compression deformation along the ½111 direction, as illustrated in the schematic (a), where averaged stresses are obtained over 2000 AIMD steps for each data point. (e,f) All six instantaneous stress components are plotted for shear and compression at 3000 K, together with two representative average (over 2000 AIMD steps) stress curves for each case. (A colour version of this figure can be viewed online.)

Fig. 1(c) and (d), we plot the calculated average (over 2000 AIMD steps) stress-strain curves of diamond under shear and compression strains at different temperatures in selected directions with the highest compression strength and the lowest shear strength, respectively, as depicted in Fig. 1(a). In our notation, we use ½hkl to ! ! ! denote a direction determined by h a 1 þ k a 2 þ l a 3 ; and use the Miller indices ðhklÞ to represent a family of crystal planes normal to ! ! ! ! ! the direction h b 1 þ k b 2 þ l b 3 , where a i and b i (i ¼ 1; 2; 3) are the direct and reciprocal lattice basis vectors, respectively. The mechanical strength, defined as the peak stress in these stressstrain curves [42e44], of diamond is reduced by about half as the temperature reaches 3000 K, indicating a much more sensitive

temperature dependence of the plastic mechanical properties compared to the elastic properties, in agreement with the experimental results [17e20]. The insensitivity of elastic constants to temperature can also be seen from the minimal differences among the slopes of the stress-strain curves near zero strain at different temperatures under shear or compression strains. Although high temperature (3000 K) causes large thermal fluctuations in instantaneous stress responses during shear and compression AIMD simulations, apart from the average shear (szx ) or compression (szz ) component, all the other five independent stresses are averaged (over 2000 AIMD steps) to zero, as shown in Fig. 1(e) and (f), as required by our AIMD simulation process described above.

364

L. Wen et al. / Carbon 155 (2019) 361e368

To better understand the temperature effect on the mechanical failure mechanism in diamond from an atomistic perspective, we examine the evolution of bond length and angle during deformation at different temperatures. As shown in the schematic in Fig. 1(a), the carbon sp3 tetrahedrons are the building blocks of the diamond structure. The straining of the supercell results in continuous distortions of the tetrahedrons across the supercell. We highlight one representative tetrahedron of the supercell (as shown in Fig. 2(a) and (d)) and calculate the bond lengths and angles. During the process of shear loading along the ð111Þ½112 slip direction, bonds OA and OD are elongated, while bonds OB and OC are compressed. The degree of distortion of bond OA is the largest among all the bonds, and we plot the length of bond jOAj against shear strain in Fig. 2(b). We also plot in Fig. 2(c) the most distorted bond angle :AOD among the six angles centered around atom O. Under compression along the ½111 direction, the distortion of bonds OA, OB and OC is symmetrical and the bond length is only slightly shortened (within 1% at all the examined temperatures), while bond OD is extensively compressed (up to 15% at 300 K before breaking). In Fig. 2(e) and (f), we plot the length of bond jODj and the most distorted bond angle :AOD. The results in Fig. 2 show that, with increasing temperature, the maximum strain at the failure point becomes smaller under both shear and compression. This temperature softening effect can be understood by noticing the temperature dependent thermal fluctuation of the bond angle :AOD, whose variation at different temperatures is shown in Fig. 2(c) and (f) for shear and compression, respectively. At higher temperatures, the thermal fluctuation of :AOD becomes larger, making it easier to overcome the kinetic barrier and cause subsequent breaking of the initial sp3 structure. In

Fig. 2(c), both the critical shear strains, at which the sp3 structure breaks, and the corresponding averaged maximum values of :AOD at these failure strains, shown by the peaks of the white curves, vary for different temperatures, but the high fluctuation probability upper bounds of the angle :AOD at the failure strains are close, at approximately 130+, indicated by the red horizontal line in Fig. 2(c). This angle is much larger than that of the initial unstrained sp3 tetrahedron (109:5+ ), and even goes beyond that of sp2 honeycomb (120+ ), resulting in an unstable tetrahedron structure and subsequent mechanical failure when large numbers of local sp3 tetrahedrons are deformed with such high fluctuation probabilities. Similar arguments can be applied to the compression case shown in Fig. 2(f), where the high fluctuation probability lower bounds of the angle :AODz92+ is close to a right-angle, which means that the four atoms OABC are nearly in the same plane, resembling an sp2 structure. This makes the fourth bond OD in the tetrahedron (see Fig. 2(d)) highly unstable, resulting in mechanical failure when large numbers of local sp3 tetrahedrons are deformed. To evaluate mechanical properties of diamond at extreme HPHT, we have calculated shear stress-strain curves at hydrostatic pressure p ¼ 300 and 600 GPa and temperature T ¼ 6500 and 7000 K. For completeness of the pressure range, we also include the results for p ¼ 100 GPa and T ¼ 5500 K. The (111)½112 shear direction we consider here is the easy-slip direction of diamond, which has been reported to have the lowest shear strength on the easy cleavage (111) plane [42], thus limiting the overall mechanical performance of diamond. In Fig. 3, we plot calculated average (over 2000 AIMD steps) stress-strain curves of diamond for these temperatures and pressures under shear loading in the ð111Þ½112 slip direction. The temperatures in our calculation are chosen to be lower than the

Fig. 2. Evolution of bond lengths and angles with shear and compression strains at different temperatures. (a,b,c) The most distorted bond length jOAj and angle :AOD, under shearing along the ð111Þ½112 slip direction. (d,e,f) The most distorted bond length jODj and angle :AOD, under compression along the ½111 direction. In all presented cases, averaged results over every 2000 AIMD steps are also plotted as curves in white color. (A colour version of this figure can be viewed online.)

L. Wen et al. / Carbon 155 (2019) 361e368

Fig. 3. Stress-strain relations of diamond under shear loading along the ð111Þ½112 slip direction at different pressures and temperatures. The dashed lines mark the shear strains at which the small-size coexistence method is employed to determine the thermodynamic stability and melting temperature of the shear strained diamond structures. (A colour version of this figure can be viewed online.)

melting points Tm of the unstrained diamond structures at the corresponding pressures by about 500 K or 1000 K, with Tm estimated using the small-size coexistence method [31,32] as depicted in Fig. 4(a) (see below). The selection of the temperatures ensures

365

that these probing regimes are near the melting line of the phase diagram of diamond, which relates the mechanical stability and the phase stability, as we will discuss later. The strength of diamond is affected by both the hydrostatic pressure p and temperature T. As can be seen in Fig. 3, the maximum shear stress is larger for higher pressure or lower temperature. Even at such high temperatures, the shear strength of diamond is still very high if large enough hydrostatic pressure exists. The results show that the shear strength, which is closely related to hardness, of diamond at extreme HPHT conditions corresponding to the core regions of Neptune and Uranus (Fig. 3) is as high as it is on Earth's surface under ambient conditions (Fig. 1(c)). It should be noted that the melting of a solid strained structure at such high temperatures may not be directly seen due to the overheating issue in the calculations, which is well known in the “heat until it melts” approach [49,50]. Previous works on strain-induced reduction of melting temperature [51e53] suggest that a strained crystal might melt at temperatures below the unstrained melting temperature Tm . We now turn to the analysis of the phase stability and shearinduced early melting of the strained diamond structure. The workflow in Fig. 4(a) shows the steps of the small-size coexistence method [31,32], which we employed to determine whether a solid state will persist or melt according to the statistical ratio a≡NLiquid =NTotal , where NTotal (¼10) is the total number of AIMD runs for the same initial solid-liquid coexisting state and NLiquid represents the number of the runs that end up in a liquid state, see Methods of Calculations above for details. Study of melting properties of diamond using the small-size coexistence method at

Fig. 4. Phase stability determined using the small-size coexistence method. (a) Diagram demonstrating the small-size coexistence method. The initial diamond structure is the one with an applied shear strain at a given temperature and pressure, taken from the results in Fig. 3. (b) MSD of carbon atoms in the solid-liquid coexistence vs. AIMD time at different temperatures under zero pressure. The atomic structures at the last AIMD steps (at t ¼ 6 ps) are shown as insets. (c) Estimation of melting temperatures of diamond with no strain (in equilibrium) under different hydrostatic pressures using the small-size approach and the bisection method in a few iteration steps. (d) Determination of the stable state (solid or liquid) of the diamond structure with applied shear strains under different temperatures and pressures. Points with a≡NLiquid =NTotal > 0:5 correspond to liquid. For all the reported calculations, we use NTotal ¼ 10. (A colour version of this figure can be viewed online.)

366

L. Wen et al. / Carbon 155 (2019) 361e368

zero pressure is complicated by the existence of a third, graphitic phase. In Fig. 4(b), the results of the (diamond and liquid carbon) coexistence tests at three temperatures are shown for the unstrained diamond at p ¼ 0 GPa. At each temperature, we perform NTotal ¼ 10 coexistence runs and choose a typical run to plot the mean squared displacement (MSD) and display the atomic structure of the last AIMD step. The MSD is defined as a function of the simulation time t,

MSDðtÞ ¼

N 1 X 2 ! ! ½ r i ðtÞ  r i ð0Þ N i¼1

(1)

! where r i ðtÞ is the position of atom i at time t and N the total ! number of atoms. Note that r i ð0Þ is the position of atom i of the initial structure of the coexistence, thus here the MSD is defined relative to the initial structure, similar to the previous definition [24]. In Fig. 4(b), the MSD at 3000 K (4000 K) increases during the first 0:5 ps (2 ps) and then saturates, where the initial increase in the MSD corresponds to the (partial) transformation of the coexistence to graphitic layers, as shown by the insets in Fig. 4(b). The saturation of the MSD suggests that the system stays in the solid (graphite or diamond) state, where the atoms are finally confined to local positions. At 4500 K, most of the coexistence runs exhibit linearly diverged MSD and the final state is liquid, thus giving an estimated range of the transition temperature Tm from solid carbon to liquid carbon at zero pressure 4000 K < Tm < 4500 K, which is consistent with the previously reported value 4160 K [54]. It should also be noted that the melting of diamond at zero pressure can only happen at the metastable continuation of the melting line, as constructed by Basharin et al. [54]. The detailed analysis of the structure melting evolution at T ¼ 4500 K shows that the diamond structure can directly melt into liquid carbon, without the formation of oriented layer structure of graphite (see the Supplementary Materials), which further confirms the consistency between our results and the previous conclusion of Basharin et al. [54]. The outcome liquid carbon has a complex composition of carbon atoms with different hybridization, which is 25%ðsp1 Þ, 65%ðsp2 Þ and 10%ðsp3 Þ, with a density of 1:99 g =cm3 . This composition generally agrees with the previous results calculated using the same method by Dozhdikov et al. for the liquid carbon with the same density at higher temperature and pressure [55]. The formation of graphitic phase at 3000  4000 K indicates that graphite is thermodynamically more stable at zero pressure in comparison with diamond, which is consistent with the fact that diamond is a metastable phase at ambient pressure. Earlier CarParrinello MD or AIMD calculations have studied graphitization of a diamond surface at zero pressure under different temperatures [56,57]. Their results show that diamond starts to transform into graphite from the top layers of its (111) surface at temperatures ranging from 2500 K to 3900 K, depending on weather the surface is reconstructed, stepped or hydrogen saturated in a vacuumdiamond slab model. However, these surface characteristics, such as reconstructions, steps and hydrogen saturations, do not exist inside bulk diamond. The graphitization temperature (3000  4000 K) we obtained for bulk diamond generally agree with those earlier studies of diamond surfaces [56,57]. But the (solid-liquid two-phase) coexistence model adopted in this work, which can be viewed as a liquid-diamond slab model, may not be applicable to studies of graphitization and shear melting of bulk diamond at zero pressure, where diamond, graphite and liquid carbon phases transform among each other. When pressure is high (p > 100 GPa), where diamond becomes the most stable solid phase, the coexistence model is expected to well predict the melting temperature of

(strained) diamond. We now focus on the high-pressure regimes. We first calculate the melting points (Tm ) of the unstrained diamond under given pressures, where we apply the workflow multiple times at a series of temperatures determined by the iteration scheme of the bisection method, until Tm is located in a narrow temperature window where a≡NLiquid =NTotal jumps from 0 to 1, indicating melting. The calculated results are plotted in Fig. 4(c), which give melting points in good agreement with those in previous works [22]. Similarly, we checked the states of the shear strained diamond under different pressures and temperatures at several chosen strains as illustrated in Fig. 3. The results are plotted in Fig. 4(d). All the unstrained (εzx ¼ 0) diamond structures at these selected pressures and corresponding temperatures have a ¼ 0 (not shown), thus being stable as solid. With increasing shear strains, the ratio a increases, some of which become larger than 0.5, meaning that these shear strained diamond solid structures are thermodynamically unstable compared to the liquid state. The estimated reduction (DTm ) of the shear-induced melting temperature relative to Tm under different pressures is: DTm  400 K for p ¼ 100 GPa; DTm  1150 K for p ¼ 300 GPa; and DTm  800 K for p ¼ 600 GPa. To understand the structure of the stable phases of shear strained diamond, along its ð111Þ½112 easy-slip direction at different pressures and temperatures, in the NpT AIMD simulation with the small-size coexistence method, two typical examples are plotted in Fig. 5 to demonstrate the stresses evolution for two independent coexistence tests with similar random initial liquid configurations but different final states at 300 GPa, 7000 K, where the applied shear strain is fixed at εzx ¼ 0:140 during the AIMD run. In the first example (Fig. 5(a)), the final state of the coexistence AIMD run with a 50% probability (see Fig. 4(d)) is in liquid state, while in the second example (Fig. 5(b)), it recovers with a 50% probability (see Fig. 4(d)) the initial solid state. The insets of the upper panels in Fig. 5 are the evolution of MSD with time, as defined in Eq. (1), indicating one liquid state and one solid state, respectively, which can also be seen from the atomic structure snapshots at the final step (t ¼ 6 ps) shown as the insets of the lower panels. After melting occurs the shear stresses (szx ) completely relax (see Fig. 5(a)), suggesting that the resulting liquid is an unstable supercooling liquid, since the given pressure and temperature is within the solid-state regime of the p -T phase diagram for unstrained diamond. The liquid may crystallize into diamond or solidify into amorphous structure [58] in real experiments, with shear stresses fully released in the newly formed solid. A similar phenomenon is previously reported by Levitas et al. for Cu and Al in a classic MD shock compression simulation [28], and termed as virtual melting. In their work, high strain rates (0:001e1ps1 ) are used to simulate the shock wave propagation and reduction of melting temperature as large as 80% has been observed for Cu and Al. We have adopted the lower limit (0:001ps1 ) of the strain rate used by Levitas et al. in our simulations. Using strain rates lower than 0:001ps1 in our AIMD simulations would have been computationally prohibitively expensive, thus not attempted in this work. Recent experiments on nanosecond pulse laser mediated phase conversion from sp2 -rich carbon into diamond also reveal similar super undercooling carbon liquid and the subsequent crystallization [41,59e61]. The virtual melting behavior caused by shear stresses reduces diamond's ideal strength in HPHT applications. For example, the extensive compression and shear deformations of the indenter surface due to high pressure in DAC [62] give rise to large shear stresses, which may greatly reduce the viable high-temperature and high-pressure regime for LHDACs. In Fig. 6, we plot our results on the carbon phase diagram to demonstrate this reduction in the melting temperature of diamond.

L. Wen et al. / Carbon 155 (2019) 361e368

367

Fig. 5. Stress evolutions for two independent coexistence tests with similar random initial liquid configurations at 300 GPa and 7000 K with an applied shear strain szx ¼ 0:140. A case with a liquid final state is shown in (a) and the other with a solid final state is shown in (b). The insets of the upper panels are the evolution of MSD with time (in unit of ps), defined in Eq. (1). The insets of the lower panels are the atomic structure snapshots at the final AIMD step (at t ¼ 6 ps). The dashed lines denote the expected average stresses (over 2000 AIMD steps), which should be i) sxx ¼ syy ¼ szz ¼ 300 GPa; ii) sxy ¼ syz ¼ 0 and iii) szx ¼ 0 for the liquid state, or 64:5 GPa for the solid state which is the shear stress of the initially shear strained structure before the coexistence test, as shown in Fig. 3. (A colour version of this figure can be viewed online.)

and analysis that plastic strengths of diamond reduce quickly as the temperature reaches above 1000 K, while its elastic constants remain almost temperature independent. At pressures higher than 100 GPa, there is an unexpectedly large reduction of melting temperature caused by a shear-induced early melting of diamond, by as much as 1150 K, under strains along the ð111Þ½112 easy-slip shear direction. This result profoundly alters the p-T phase diagram and the equation of state of diamond exposing to shear stresses, with significant implications for LHDAC or shock wave compression experiments where large shear stresses commonly exist. These findings establish new benchmarks and provide crucial insights for understanding the extreme mechanics of diamond under HPHT loading conditions. These benchmarks and insights offer important new ingredients for consideration in the design and operation of experiments employing DAC devices under extreme loadings, and also introduce new constraints for modeling the behavior of diamond as a constituent component in carbon-rich giant planetary interiors. Fig. 6. Shift of the melting line of diamond on the phase diagram of carbon. The black dashed lines are the phase boundaries of carbon under hydrostatic pressures determined by Correa et al. [22], which is consistent with our results represented by the red circles. The thick line denotes the estimated upper bound of the melting temperature for large-shear-strained diamond crystal along the shear direction ð111Þ½112, calculated under different shear strains indicated by different colored triangles. (A colour version of this figure can be viewed online.)

As can be seen from the phase diagram, a significant fraction of the solid diamond regime is predicted to be unstable, if large shear strain exists.

Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant No. 11574197) and Ministry of Science and Technology of China (Grant No. 2016YFA0300500). Computation was performed at the Center for High Performance Computing, Shanghai Jiao Tong University.

Appendix A. Supplementary data 4. Conclusions In summary, we have shown by systematic AIMD simulations

Supplementary data to this article can be found online at https://doi.org/10.1016/j.carbon.2019.08.079.

368

L. Wen et al. / Carbon 155 (2019) 361e368

References [1] A. Salamat, R.A. Fischer, R. Briggs, M.I. McMahon, S. Petitgirard, In situ synchrotron X-ray diffraction in the laser-heated diamond anvil cell: melting phenomena and synthesis of new materials, Coord. Chem. Rev. 277e278 (2014) 15e30. [2] A. Zerr, et al., Synthesis of cubic silicon nitride, Nature 400 (1999), 340-340. [3] E. Gregoryanz, et al., Synthesis and characterization of a binary noble metal nitride, Nat. Mater. 3 (2004), 294-294. [4] M.I. Eremets, A.G. Gavriliuk, I.A. Trojan, D.A. Dzivenko, R. Boehler, Singlebonded cubic form of nitrogen, Nat. Mater. 3 (2004), 558-558. [5] V.L. Solozhenko, O.O. Kurakevych, D. Andrault, Y. Le Godec, M. Mezouar, Ultimate metastable solubility of boron in diamond: synthesis of superhard diamondlike BC5, Phys. Rev. Lett. 102 (2009), 015506-015506. [6] L. Zhang, et al., Disproportionation of (Mg,Fe)SiO3 perovskite in Earth's deep lower mantle, Science 344 (2014), 877-877. [7] K. Hirose, et al., Crystallization of silicon dioxide and compositional evolution of the Earth's core, Nature 543 (2017), 99-99. [8] K. Ohta, Y. Kuwayama, K. Hirose, K. Shimizu, Y. Ohishi, Experimental determination of the electrical resistivity of iron at Earth's core conditions, Nature 534 (2016), 95-95. ^pkova , R.S. McWilliams, N. Go mez-Pe rez, A.F. Goncharov, Direct [9] Z. Kono measurement of thermal conductivity in solid iron at planetary core conditions, Nature 534 (2016), 99-99. [10] E. Anders, E. Zinner, Interstellar grains in primitive meteorites: diamond, silicon carbide, and graphite, Meteoritics 28 (1993) 490e514. [11] S. Richter, U. Ott, F. Begemann, Tellurium in pre-solar diamonds as an indicator for rapid separation of supernova ejecta, Nature 391 (1998) 261e263. [12] Z.R. Dai, et al., Possible in situ formation of meteoritic nanodiamonds in the early Solar System, Nature 418 (2002) 157e159. [13] M. Ross, The ice layer in Uranus and Neptune–diamonds in the sky? Nature 292 (1981), 435-435. [14] S. Brygoo, et al., Laser-shock compression of diamond and evidence of a negative-slope melting curve, Nat. Mater. 6 (2007), 274-274. [15] J.H. Eggert, et al., Melting temperature of diamond at ultrahigh pressure, Nat. Phys. 6 (2009), 40-40. [16] R.F. Smith, et al., Ramp compression of diamond to five terapascals, Nature 511 (2014), 330-330. [17] E.S. Zouboulis, M. Grimsditch, A.K. Ramdas, S. Rodriguez, Temperature dependence of the elastic moduli of diamond: a Brillouin-scattering study, Phys. Rev. B 57 (1998) 2889e2896. [18] D.R. Mumm, K.T. Faber, M.D. Drory, C.F. Gardinier, High-temperature hardness of chemically vapor-deposited diamond, J. Am. Ceram. Soc. 76 (1993) 238e240. [19] N.V. Novikov, Y.V. Sirota, V.I. Mal'nev, I.A. Petrusha, Mechanical properties of diamond and cubic BN at different temperatures and deformation rates, Diam. Relat. Mater. 2 (1993) 1253e1256. [20] V.A. Mukhanov, O.O. Kurakevych, V.L. Solozhenko, Hardness of materials at high temperature and high pressure, Philos. Mag. 89 (2009) 2117e2127. [21] X. Wang, S. Scandolo, R. Car, Carbon phase diagram from ab initio molecular dynamics, Phys. Rev. Lett. 95 (2005), 185701-185701. [22] A.A. Correa, S.A. Bonev, G. Galli, Carbon under extreme conditions: phase boundaries and electronic properties from first-principles theory, Proc. Natl. Acad. Sci. U.S.A. 103 (2006), 1204-1204. [23] M. Martinez-Canales, C.J. Pickard, R.J. Needs, Thermodynamically stable phases of carbon at multiterapascal pressures, Phys. Rev. Lett. 108 (2012), 045704-045704. [24] C. Eisenmann, C. Kim, J. Mattsson, D.A. Weitz, Shear melting of a colloidal glass, Phys. Rev. Lett. 104 (2010), 035502-035502. [25] C.J. Wu, P. Soderlind, J.N. Glosli, J.E. Klepeis, Shear-induced anisotropic plastic flow from body-centred-cubic tantalum before melting, Nat. Mater. 8 (2009) 223e228. [26] C. Healy, S. Koch, C. Siemers, D. Mukherji, G.J. Ackland, Shear melting and high temperature embrittlement: theory and application to machining titanium, Phys. Rev. Lett. 114 (2015), 165501-165501. [27] Y. Feng, J. Goree, B. Liu, Evolution of shear-induced melting in a dusty plasma, Phys. Rev. Lett. 104 (2010), 165003-165003. [28] V.I. Levitas, R. Ravelo, Virtual melting as a new mechanism of stress relaxation under high strain rate loading, Proc. Natl. Acad. Sci. U.S.A. 109 (2012) 13204e13207. [29] M.M. Budzevich, V.V. Zhakhovsky, C.T. White, I.I. Oleynik, Evolution of shockinduced orientation-dependent metastable states in crystalline aluminum, Phys. Rev. Lett. 109 (2012), 125505-125505. [30] S. Butler, P. Harrowell, Factors determining crystal-liquid coexistence under shear, Nature 415 (2002), 1008-1008. [31] Q.-J. Hong, A. van de Walle, Solid-liquid coexistence in small systems: a statistical method to calculate melting temperatures, J. Chem. Phys. 139 (2013), 094114-094114. [32] Q.-J. Hong, A. van de Walle, A user guide for SLUSCHI: solid and liquid in ultra

small coexistence with hovering interfaces, Calphad 52 (2016) 88e97. [33] G. Kresse, Hafner, J. Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47 (1993) 558e561. [34] G. Kresse, Hafner, J. Ab initio molecular-dynamics simulation of the liquidmetal–amorphous-semiconductor transition in germanium, Phys. Rev. B 49 (1994) 14251e14269. [35] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (1996) 11169e11186. [36] G. Kresse, D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59 (1999) 1758e1775. [37] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77 (1996) 3865e3868. [38] H.J. Monkhorst, J.D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B 13 (1976) 5188e5192. [39] M. Parrinello, A. Rahman, Crystal structure and pair potentials: a moleculardynamics study, Phys. Rev. Lett. 45 (1980) 1196e1199. [40] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 1987. [41] J. Narayan, A. Bhaumik, S. Gupta, A. Haque, R. Sachan, Progress in Q-carbon and related materials with extraordinary properties, Mater. Res. Lett. 6 (2018) 353e364. [42] Y. Zhang, H. Sun, C. Chen, Atomistic deformation modes in strong covalent solids, Phys. Rev. Lett. 94 (2005), 145505-145505. [43] Z. Pan, H. Sun, Y. Zhang, C. Chen, Harder than diamond: superior indentation strength of wurtzite bn and lonsdaleite, Phys. Rev. Lett. 102 (2009), 055503055503. [44] B. Li, H. Sun, C. Chen, Extreme mechanics of probing the ultimate strength of nanotwinned diamond, Phys. Rev. Lett. 117 (2016), 116103-116103. [45] Q.-J. Hong, S.V. Ushakov, A. Navrotsky, A. van de Walle, Combined computational and experimental investigation of the refractory properties of La2Zr2O7, Acta Mater. 84 (2015) 275e282. [46] Q.-J. Hong, et al., Combined computational and experimental investigation of high temperature thermodynamics and structure of cubic ZrO2 and HfO2, Sci. Rep. 8 (2018), 14962-14962. [47] J. Xie, S.P. Chen, J.S. Tse, S.d. Gironcoli, S. Baroni, High-pressure thermal expansion, bulk modulus, and phonon structure of diamond, Phys. Rev. B 60 (1999) 9444e9449. ~ ez Valdez, K. Umemoto, R.M. Wentzcovitch, Elasticity of diamond at [48] M. Nún high pressures and temperatures, Appl. Phys. Lett. 101 (2012), 171902171902. [49] J.-Y. Raty, E. Schwegler, S.A. Bonev, Electronic and structural transitions in dense liquid sodium, Nature 449 (2007), 448-448. rah, Melting curve of aluminum up to 300 [50] J. Bouchet, F. Bottin, G. Jomard, G. Ze GPa obtained through ab initio molecular dynamics simulations, Phys. Rev. B 80 (2009), 094102-094102. [51] U. Tartaglino, E. Tosatti, Strain effects at solid surfaces near the melting point, Surf. Sci. 532e535 (2003) 623e627. [52] T. Frolov, Y. Mishin, Effect of nonhydrostatic stresses on solid-fluid equilibrium. I. Bulk thermodynamics, Phys. Rev. B 82 (2010), 174113-174113. [53] Y.S. Hwang, V.I. Levitas, Internal stress-induced melting below melting temperature at high-rate laser heating, Appl. Phys. Lett. 104 (2014), 263106263106. [54] A.Y. Basharin, V.S. Dozhdikov, A.V. Kirillin, M.A. Turchaninov, L.R. Fokin, Phase diagram with a region of liquid carbon-diamond metastable states, Tech. Phys. Lett. 36 (2010) 559e562. [55] V.S. Dozhdikov, A.Y. Basharin, P.R. Levashov, D.V. Minakov, Atomistic simulations of the equation of state and hybridization of liquid carbon at a temperature of 6000 K in the pressure range of 1-25 GPa, J. Chem. Phys. 147 (2017), 214302-214302. [56] A. De Vita, G. Galli, A. Canning, R. Car, A microscopic model for surfaceinduced diamond-to-graphite transitions, Nature 379 (1996) 523e526. [57] G. Kern, J. Hafner, Ab initio molecular-dynamics studies of the graphitization of flat and stepped diamond (111) surfaces, Phys. Rev. B 58 (1998) 13167e13175. [58] G. Moras, et al., Shear melting of silicon and diamond and the disappearance of the polyamorphic transition under shear, Phys. Rev. Mater. 2 (2018), 083601-083601. [59] J. Narayan, A. Bhaumik, Novel phase of carbon, ferromagnetism, and conversion into diamond, J. Appl. Phys. 118 (2015), 215303-215303. [60] J. Narayan, et al., Direct conversion of carbon nanofibers and nanotubes into diamond nanofibers and the subsequent growth of large-sized diamonds, Nanoscale 11 (2019) 2238e2248. [61] A. Haque, R. Sachan, J. Narayan, Synthesis of diamond nanostructures from carbon nanotube and formation of diamond-CNT hybrid structures, Carbon 150 (2019) 388e395. [62] B. Li, et al., Diamond anvil cell behavior up to 4 Mbar, Proc. Natl. Acad. Sci. U.S.A. 115 (2018) 1713e1717.