Tensorial stress−strain fields and large elastoplasticity as well as friction in diamond anvil cell up to 400 GPa

Various phenomena (fracture, phase transformations, and chemical reactions) studied under extreme pressures in diamond anvil cell are strongly affected by fields of all components of stress and plastic strain tensors. However, they could not be measured. Here, we suggest a coupled experimental−theoretical−computational approach that allowed us (using published experimental data) to refine, calibrate, and verify models for elastoplastic behavior and contact friction for tungsten (W) and diamond up to 400 GPa and reconstruct fields of all components of stress and large plastic strain tensors in W and diamond. Despite the generally accepted strain-induced anisotropy, strain hardening, and path-dependent plasticity, here we showed that W after large plastic strains behaves as isotropic and perfectly plastic with path-independent surface of perfect plasticity. Moreover, scale-independence of elastoplastic properties is found even for such large field gradients. Obtained results open opportunities for quantitative extreme stress science and reaching record high pressures.


INTRODUCTION
In static high-pressure research, megabar pressures are generated by compression of a thin sample by two diamonds in diamond anvil cells (DAC) [1][2][3][4] ; see Fig. 1. This process is accompanied by large plastic deformation of a sample and large elastic deformation of the diamond. 5,6 Various problems, such as the study of physical, chemical, geological, and mechanical phenomena and synthesis of new phases in a sample, as well as the increasing range of achievable pressures,  are related to knowledge of the fields of all components of the stress, elastic, and plastic strain tensors in DAC. While most measurements and discussions are related to pressure only, it is evident that elastic deformation and fracture of diamond and plastic flow of a sample and gasket depend on all components of the stress tensor. Contact friction between diamond and sample/gasket plays a key role in generating high pressure without fracture of the diamond; friction is a shear stress that depends on the stress normal to the contact surface. It is also well-known that phase transformations and chemical reactions in solids depend not only on pressure, but also on the deviatoric stresses and plastic strains. 12,[14][15][16][20][21][22][23] All of these fields are extremely complex and heterogeneous, e.g. with normal stresses varying by megabar over 20 μm. 6,7 Measurement of the radial pressure distribution at the samplediamond boundary was based on the ruby fluorescence method, which worked up to 185 GPa. 4 For higher pressure, radial pressure distribution averaged over the sample thickness is determined using X-ray diffraction in a sample. 6,7 The radial thickness profile, which characterizes both elastic deformation of an anvil and elastoplastic deformation of a sample/gasket, was measured utilizing in situ high-pressure X-ray absorption. 6,7 Measurement of the deviatoric stress was limited to the difference between axial stresses σ zz and radial stresses σ rr averaged over the entire sample. 8,9,[17][18][19] Plastic deformation fields in the sample compressed in DAC and contact friction stresses were not measured at all. Thus, despite significant progress, it is unlikely that all tensorial fields in DAC will be measured. Theoretical approaches and finite element method (FEM) simulations [23][24][25][26][27][28][29] of the DAC are based on relatively simple models with linear pressure dependence of the yield strength and simplified contact friction conditions. The most sophisticated model and the best numerical reproduction of the experimental pressure distribution in ref. 5 was obtained in ref. 27 for compression of rhenium up to 285 GPa. However, in that study the plastic sliding along the contact surface and also the dependence of the friction coefficient on the normal contact pressure were ignored. Also, good correspondence was obtained for one pressure distribution only; for two smaller pressure levels, significant deviation from the experiment existed, i.e. the model is not adequate. To obtain such a description of the experiment, the third-order elastic constants of diamond were modified. Besides, since experimental thickness profiles of the sample were not available, reproducing the pressure distributions was aimed only. Thus, while essential improvement in reproducing pressure distribution in comparison with the previous works 26,27 was achieved, obtained mechanical properties cannot be considered as verified, and corresponding stress and plastic strain tensor fields may contain significant inaccuracies.
We suggest the following coupled experimental−theoretical −computational approach for determination of all stress and plastic strain tensorial fields, elastoplastic properties, and contact friction rules. All fields that can be measured should be measured. Physics-based models for elastoplastic behavior and contact friction should be iteratively developed and refined, and all material properties should be calibrated by fitting to some experimental fields and verified by comparison with other experimental fields. With these properties, simulations provide all fields, including components of the stress and plastic strain tensors, friction stress, etc., i.e., those which cannot be directly measured. To obtain the first results from this method, we will use the most advanced experimental data on compression of W in DAC up to 400 GPa 6 and slightly generalize our models for large elastoplastic deformations and contact friction from. 24,27 RESULTS AND DISCUSSION Model A complete system of equations for fourth-order elasticity of diamond, large elastoplastic deformation of W, combined Coulomb and plastic friction, geometry of DAC, formulation of axisymmetric problem in cylindrical coordinates rzθ, and nonlinear elastic properties are presented in the Methods section. It is known that the yield surface in the six-dimensional space of components of the stress tensor evolves during plastic deformation ( Fig. 13), exhibiting strain hardening; this evolution depends on the entire history of plastic strain, and material acquired deformation-induced anisotropy. 24,30 It was suggested in ref. 24 as the postulate of the perfect plasticity that, above some level of accumulated plastic strain q > m (q is defined in Eq. (11)) and for a deformation path without sharp changes in directions (monotonous deformation), the initially isotropic polycrystalline materials are deformed as perfectly plastic and isotropic with a strainhistory-independent surface of the perfect plasticity (Fig. 13). This statement means that (1) the strain hardening is saturated, and that (2) strain-induced anisotropy and path dependence do not exhibit themselves at monotonous loading. Some qualitative arguments in favor of the postulate of the perfect plasticity have been analyzed, 24 but quantitative proof was not given for any material. Here, we incorporated this postulate into our model and will prove that such a model describes well experimental data. Our model is based on the linear pressure dependence of the yield strength in compression σ y ¼ σ 0y þ ap with two material parameters, with no plastic strain or plastic strain path dependence. Another hypothesis that will be proven is that despite the μmsized sample thickness and huge stress and plastic strain gradients, i.e. conditions that require utilization of scaledependent and the gradient plasticity, 31-33 much simpler local plasticity provides adequate description of experiments.
Contact friction stress is determined either by the Coulomb law where σ c is the normal contact stress and μ is the friction coefficient, or by the yield strength in shear τ f ¼ τ y ðpÞ ¼ σ y ðpÞ= ffiffi ffi 3 p (plastic friction), whichever is smaller. Sticking occurs if the contact shear stress is smaller than these critical values τ f . All of our assumptions for the yield strength (independence of plastic strain and its path) are also involved in the assumption for plastic friction. The friction coefficient is usually taken as a constant 27 because no experimental data under high pressure are available. We assume μ ¼ μ 0 þ cσ c with two material parameters.
In addition, some of the third-order elastic constants of W and fourth-order elastic constants of diamond, which are not well defined from the literature, are refined by comparison with DAC pressure and sample thickness distributions.
To summarize, in comparison with ref. 27 , current model includes fourth-order elasticity of diamond, combined Coulomb and plastic contact sliding, and linear pressure dependence of the Coulomb friction coefficient. Moreover, all unknown material parameters are calibrated using one set of experimental data and verified using another experimental set.

Results
All four material parameters in the pressure dependence of the yield strength and friction coefficient were calibrated by minimizing the error between experimental and FEM results for pressure distributions for two curves with maximum pressures p max = 170 and 240 GPa (Fig. 2a). This led to σ y ¼ 1:8 þ 0:1p; p 225 GPa μ ¼ 0:05 þ 0:001σ c ; σ c 37 GPa: (1) Unexpected strong limitations on pressure and contact stress appear because we found in FEM solutions that Coulomb sliding and plastic flow do not occur for σ c > 37 GPa and p > 225 GPa, respectively. With properties in Eq. (1), good correspondence is obtained for two other pressure distributions with p max = 300 and 400 GPa, with a maximum difference not exceeding 10% (Fig. 2a). In addition, the profile of the sample after very large compression and deformed anvil surface were reproduced for all four pressures, with a maximum discrepancy smaller than 1 μm (Fig. 2b, c). Both discrepancies are within error for an experiment under such extreme conditions. As seen in Fig. 2, properties in Eq. (1) result in having good agreement with experimental pressure distribution not only at large pressures, but also at low pressures where the error in the experimental results is assumed to be the least. The achieved good agreement at lower pressures is missing in ref. 27 . Besides, thickness profiles are properly reproduced. The curves in Fig. 2 are nontrivial, and coincidence demonstrates strong verification of the entire model and the specific material properties from Eq. (1). It also proves the validity of the postulate of the perfect plasticity for W, which was directly incorporated in our model, and sufficiency of the local elastoplastic model even at micron scale and with huge stress and plastic strain gradients. In summary, elastoplasticity and, consequently, plastic friction under such large strain and pressure is plastic strain-, plastic strain path-, and scale-independent, which drastically simplifies theory and measurements.
In addition, the higher-order elastic constants of W and diamond, which have large scatter in literature (see the Methods section), have been also refined/identified. Thus, we found the third-order constants for W, m = −1081 and n = −1164 GPa, to obtain a slightly better fit to the experimental pressure distribution curves for three lowest pressure. The fourth-order elastic constant of diamond, C 1112 = 31,214, C 1122 = 20,044, and C 1266 = 819 GPa, were found from the best fit to the sample profile under highest pressure only under constraint that they satisfy the known equation of state of diamond; see the Methods section.
The suggested method has high throughput features, which allows to determine ten material parameters using three pressure and one sample thickness distributions.
Known 5 pressure dependence of the yield strength for W has huge scatter (Fig. 3), which is related to numerous assumptions for the determination of σ y (p) and to attribution of the dependence of σ y on plastic strain to the pressure dependency. In our curve, the effect of plastic strain is excluded and the correctness of Eq. (1) is confirmed by numerous data in Fig. 2. After proving its validity, the model is used for computational reconstruction of all fields of interest. Distribution of shear friction stress and normalized radial sample velocity along the diamondsample contact surface at different pressures is shown in Fig. 4a, b. Such a complex profile of shear stresses and their evolution are nontrivial and counterintuitive. In particular, shear stress in the sticking zone makes several oscillations in a central cup region, and the sticking zone grows with increasing compression. The plastic friction zone is surprisingly narrow, which does not allow use of the traditional method for determination of τ y ðpÞ based on a pressure gradient. 7,13,24 The maximum yield strength in shear and corresponding p in the plastic sliding zone reduces from 5.85 GPa and 77.2 GPa for p max = 164 GPa to 3.7 GPa and 44 GPa for p max = 380 GPa. The maximum shear stress in the Coulomb sliding zone is 3.21 GPa, corresponding to σ c = 37 GPa; for p max = 380 GPa, it is 2 GPa, corresponding to σ c = 26.2 GPa. An important conclusion is that, due to significant increase in the sticking zone, an increase in p max does not lead to an increase in the maximum range of σ c and friction stress, either for Coulomb or plastic friction. The only way to increase these ranges is to use torsion under a fixed force in rotational DAC, 12,15,22,[34][35][36] for which FEM simulations 29,37 show that the sticking zone is localized near the center.
The sample particles' radial velocity along the diamond-sample contact surface (Fig. 4b) is directed toward the center in the sticking zone for any pressure, and is equal, by the definition of sticking, to the velocity of the diamond contact particles. Outside the sticking zone, sample particles move away from the center, achieving maximum velocity at the edge of the culet. The maximum velocity increases to p max = 231 GPa, then reduces due to the increasing sticking zone.
All relevant fields in the central part of the W sample at maximum pressure of 300 GPa are presented in Fig. 5 on a quarter of the sample, due to the symmetries. While axial stress σ zz is independent of the z coordinate, radial stress σ rr visibly depends on z and pressure p is, by definition, in between.
All components of plastic strain and q are very heterogeneous and reach very large values. Plastic shear strain ε rz p (defined in Eq. (6)) changes sign three times in the central zone. Accumulated plastic strain q reaches its maximum value at the contact surface, especially where the thickness is smallest. Note that, for uniaxial compression/tension, q reduces to the logarithmic strain, and maximum q = 5.77 in Fig. 5c corresponds to the ratio of the initialto-final length of exp(5.77) = 321. With increasing radius, q increases further. Material rotation in Fig. 5c, which leads to the development of texture, is also very large, with a maximum of 46.8°in this region. Thus, if strain-induced anisotropy would be  (1) obtained using experimental data from ref. 6 ; symbols are from ref. 5 Large scatterin data from ref. 5 is related to multiple assumptions for the determination of the yield strength and to attribution of the dependence of σ y on plastic strain to the pressure dependency V.I. Levitas et al. present, isotropic flow theory would not describe experiments. The rotation angle, similar to shear stress σ rz , is zero at the symmetry axis and plane and increases with increasing r and z. Radial velocity (Fig. 5d) at such a pressure is directed toward the center in the entire region. It is independent of z and its magnitude increases with r. The rate of accumulated plastic strain _ q also increases with r, with zero region to the left of the white line in Fig. 5d, where plastic flow stops and the material deforms elastically. Evolution of the elastic zone with increasing pressure is shown in Fig. 2a. It appears at p max = 200 GPa and increases with increasing pressure due to cupping of diamond.
Similar fields for three other maximum pressures, which we used for calibration and verification of the model, are given in Figs 6-8. For lower pressures, we present components of the stress tensor and pressure only in Fig. 9.
In the maximum pressure range from 170 to 380 GPa, axial stress σ zz is independent of the z coordinate, radial stress σ rr shows visible dependence on z and pressure p is in between. Degree of heterogeneity of σ rr and p increases with reducing maximum pressure because of more intense plastic flow. With further maximum pressure reduction, σ zz acquires some heterogeneity at the center of the sample and around corner point between beveled and initially flat central part. The heterogeneity of σ zz is quite pronounced in the entire central region for maximum pressure in the range of 1 to 10 GPa, while below the beveled part of the anvil σ zz and even pressure is quite homogeneous along z coordinate. Shear stresses have generally similar patterns for any maximum pressure, governed by zero stress at the symmetry plane and axis and increasing shear stress with increasing z and r. In the pressure range from 1 to 10 GPa, maximum shear stress is smaller than the yield strength in shear and either sticking or Coulomb friction are involved. Due to increase in the yield strength, maximum shear stress increases when maximum pressure growth from 10 to 100 GPa, then it reduces due to the  Fields of plastic strain tensor ε p , accumulated plastic strain q, and particles' rotation angle θ, as well as normalized radial velocity v rr and _ q have similar patterns for all maximum pressures in the range from 170 to 380 GPa. However, region with elastic deformations without plasticity at the center of the sample is getting visible at p max = 240 GPa and increases with further loading.
All stress fields in the central part of the diamond for p max = 300 GPa are presented in Fig. 10. All normal stresses have their maximum at the center of the culet, with σ max zz ¼ À321 GPa and σ max rr ¼ σ max θθ ¼ À260 GPa, i.e. nonhydrostaticity is very high. Maximum shear stress σ max rz ¼ 37:5 GPa is located away from the culet. This value is significantly smaller than the theoretical shear strength of 96.6 GPa at zero pressure, which grows with pressure. 36 It is important that the regions in which maximum normal and shear stresses occur do not overlap.
The obtained fields of all components of the stress tensor are the basis for the development of criteria for fracture of diamond. To illustrate the concept, consider experimentally observed fracture due to compression stress σ [110] Fig. 10, is then σ ½110 eq ¼ σ ½110 =σ ½110 th , and fracture occurs at Fig. 6 Stress, plastic strain, and rate fields for p max = 376 GPa. a Fields of components of the stress tensor and pressure p, b plastic strain tensor ε p , c accumulated plastic strain q, and particles' rotation angle θ, and d normalized radial velocity v rr and _ q in the central part of a sample for r < 60 μm Fig. 7 Stress, plastic strain, and rate fields for p max = 240 GPa. a Fields of components of the stress tensor and pressure p, b plastic strain tensor ε p , c accumulated plastic strain q, and particles' rotation angle θ, and d normalized radial velocity v rr and _ q in the central part of a sample for r < 60 μm V.I. Levitas et al. σ ½110 eq ¼ 1. Since maximum σ ½110 eq ¼ 0:32, there is still a significant safety factor for ideal diamond along the [110] direction. For complete fracture analysis, similar distributions should be obtained for other possible fracture planes and shear stresses along these planes should be also taken into account. This is the key problem, the solution of which will allow optimization of the design of anvils and loading conditions for a perfect crystal, which will provide the upper bound of achievable pressure (stresses), provided that plastic deformation in diamond does not occur. Introducing defects into simulations will open the possibility of developing the fracture mechanics and plasticity of real diamond crystals under extreme stresses.
In summary, we suggested a novel coupled experimental −theoretical−computational approach that allowed us (using known experimental data from ref. 6 ) to extract complete information about elastoplastic properties and friction rules, as well as all complex tensorial fields for materials compressed in a DAC under extreme pressure. In particular, we refined, calibrated, and verified models for elastoplastic behavior of a sample and contact friction for W up to 400 GPa and reconstruct fields of all components of stress and large plastic strain tensors in W and diamond. In addition to quantitative information on the pressure dependence of the yield strength and friction, as well as higherorder elastic constants, we justify some general unique properties of elastoplastic behavior under very large strains and pressures: (a) Despite the generally accepted strain-induced anisotropy, strain hardening, and path-dependent plasticity, W after large plastic strains behaves isotropically and does not exhibit strain hardening and path dependence. (b) Despite the μm-sized sample thickness and huge stress (5 GPa/μm) and plastic strain gradients, scale-independence of elastoplastic properties is found.
Both of these properties drastically simplify plasticity theory and measurements under extreme conditions.
Our finding for plasticity also implies important properties for plastic friction under such extreme loading: Plastic friction is plastic strain-, plastic strain path-, and scale-independent.
Complex and counterintuitive profile of contact shear stresses and their evolution are found. The plastic friction zone is surprisingly narrow, which does not allow use of the traditional method for determination of τ y ðpÞ based on a pressure gradient. 7,13,24 Despite the maximum pressure of 380 GPa, plasticity and plastic sliding occur below 225 GPa, and Coulomb friction takes place below 37 GPa only. Due to significant increase in the sticking zone, an increase in p max above 240-300 GPa leads to decrease in the maximum range of pressure and contact stress σ c at which plastic flow, Coulomb or plastic friction occur. The only way to increase these ranges for characterization of plastic flow and contact friction is to use torsion under a fixed force in rotational DAC, 12,15,22,[34][35][36] for which FEM simulations 29,37 show that the sticking zone is localized near the center.
The field of all components of the stress tensor in diamond are the basis for the development of criteria for fracture of diamond. We illustrated the concept by considering fracture due to compression along one of the experimentally observed directions. This is an important step which after detailed study of fracture of anvil will allow optimization of the design of anvils and loading conditions for further increase in achievable pressure.
The developed experimental−theoretical−computational approach allows different realizations for different available experimental data. In particular, recently developed nanoscale sensing platform, 40 which integrates nitrogen-vacancy color centers directly into the culet of diamond anvils, allows experimental determination of distribution of all six components of the stress tensor in diamond along the flat culet, i.e., at the contact surface with the sample. This information allowed us to reconstruct nontrivial normal and shear contact stresses between diamond and gasket and fields of all components of stress tensor in the entire anvil. Since linear response is utilized, this method was applied up to 30 GPa. Fig. 8 Stress, plastic strain, and rate fields for p max = 170 GPa. a Fields of components of the stress tensor and pressure p, b plastic strain tensor ε p , c accumulated plastic strain q, and particles' rotation angle θ, and d normalized radial velocity v rr and _ q in the central part of a sample for r < 60 μm V.I. Levitas et al. Note that W is used as a gasket material in DAC at megabar pressures, i.e. obtained results have also applied importance for study of various sample materials within W gasket. Knowledge of the distributions of all (generally 12) components of stress and plastic strain tensors in a sample will allow study of their (instead of pressure alone) effect on phase transformations, chemical reactions, 12,[14][15][16][20][21][22][34][35][36][37] and various physical properties. In comparison with research under hydrostatic pressure, this will add up to 11 new dimensions to the parametric space for studying these processes, searching for new phases and materials, drastically reducing the required pressure for synthesis of new and known materials with unique properties, and understanding processes in the deep interiors of the Earth and other planets. Obtained results will also enable calibration and verification of known and new methods for measurement of the components of stress tensors in anvils and samples.

METHODS
A complete system of equations for fourth-order elasticity of diamond, large elastoplastic deformation of W, and combined Coulomb and plastic

Geometry and boundary conditions
Axisymmetric problem formulation is considered. Due to symmetry of the Mao-type DACs used in ref. 1 , only the upper part of the DAC and sample will be used in simulations. Geometry of the sample and the anvil, as well as the boundary conditions, are shown in Fig. 11 and are as follows: (1) A uniform vertical displacement is applied at the boundary between the top inclined surface of the anvil and Boehler-type seat (line CD). Distribution of stresses or displacements along this surface does not affect fields close to the diamond culet (line AG). (2) At the symmetry axis r = 0 (line AB), shear stress σ rz and horizontal displacements are zero. At the symmetry plane z = 0, shear stress σ rz and vertical displacement are zero. (3) At the contact surface between the gasket and the anvil, a combined Coulomb friction and plastic friction model, which is described below, is utilized. (4) Other surfaces not mentioned above are stress-free. Finite element algorithms for solution of the boundary-value problems are presented in Feng et al. 28 Friction model According to the combined Coulomb friction and plastic friction model, there is complete cohesion between the contact pairs unless the shear (friction) stress reaches the critical value: When friction stress reaches τ crit , contact sliding occurs in the radial direction. The critical shear stress τ crit ¼ μ σ c ð Þσ c is related to the Coulomb friction, where μ is the friction coefficient and σ c is the normal contact stress. However, the Coulomb friction stress cannot exceed the yield strength in shear τ y p ð Þ, which is defined in terms of the yield strength under compression, σ y , by τ y ¼ σ y = ffiffi ffi 3 p , based on the von Mises yield criterion. Thus, plastic sliding occurs when the Coulomb friction exceeds τ y p ð Þ. In fact, it represents plastic shear flow within a very thin material layer immediately bellow the contact surface.
In this study the yield strength and the friction coefficient are assumed to be pressure and contact pressure dependent, respectively.
We assume μ ¼ μ 0 þ cσ c with two material parameters, which after calibration, looks like μ ¼ 0:05 þ 0:001σ c ; σ c 37 GPa: Limitations on the contact stress exist because in FEM solutions, Coulomb sliding does not occur for σ c > 37 GPa, even for the highest maximum pressure of 380 GPa.

Elastoplastic material model under large strains and high pressure
We designate single and double contractions of the second-order tensors A¼ A ij È É and B¼ B ij È É over one and two indices as A Á B¼ A ij B jk È É and A : B¼ A ij B ji È É , respectively. The subscript 's' denotes symmetrization, and the subscripts 'e' and 'p' denote elastic and plastic part of a tensor, respectively. The superscripts −1 and 'T' designate the inverse and transposition of a tensor. I is the second-order unit tensor.
The complete system of equations for a large elastoplastic deformation of a sample is as follows: 25,28 Decomposition of the deformation gradient F in to elastic F e and plastic F p parts where r and r 0 are the position vectors of material points in the actual (deformed) configuration and the reference (undeformed) configuration, respectively; V e and V p are symmetric elastic and plastic left stretch tensors, respectively, U p is the plastic right stretch tensor, and R e is the proper orthogonal elastic rotation tensor. Elastic strain B e and its Jaumann objective time derivative Fig. 10 Field of the components of the stress tensor, pressure, and equivalent stress σ ½110 eq near the tip of a diamond anvil for r < 100 and z < 70 μm for p max = 300 GPa Decomposition of the velocity gradient l, into symmetric deformation rate d and skew symmetric spin w where D p is the plastic deformation rate. Isotropic elasticity rule Here σ is the true Cauchy stress, J = def F is the Jacobian, and Ψ is the specific Helmholtz free energy per unit undeformed volume.
Pressure-dependent yield surface (surface of perfect plasticity) where s is the deviatoric part of Cauchy stress σ, and σ y is the yield strength in compression. Plastic flow rule where λ ≥ 0 is a scalar determined from the consistency condition _ φ ¼ 0. The rate of accumulated plastic strain (plotted in Figs 5-8) Equilibrium condition where ∇· is the divergence operator in the deformed configuration.

The yield strength in compression
We assume σ y ¼ σ 0y þ ap with two material parameters, which, after calibration, results in σ y ¼ 1:8 þ 0:1p; p 225 GPa: Limitation on the pressure exists because, in FEM solutions, plastic flow does not occur for p > 225 GPa, despite the maximum pressure of 380 GPa.

Material properties
Diamond. All elastic material constants are taken from Telichko et al., 42 which, to the authors' knowledge, is the only reference that provides all third-and fourth-order elastic constants for diamond. These were determined using first principle simulations. Thus, we used the following elastic constants in our simulations: Since available data for the third-order elastic constants from different references have significant scatter, [43][44][45][46] we assume that some of the fourth-order elastic constants are not precise either. Indeed, for the elastic constants from Telichko et al., 42 we were unable to obtain the experimental equation of state collected in Maezono et al. 47 . Thus, we changed C 1112 , C 1122 , and C 1266 to the values indicated in Eq. (21) in order to receive good correspondence with the equation of state from Sato et al. 48 ; see Fig. 12, and sample profile at highest pressure, see Fig. 2c.
A majority of equations of state of diamond determined by different methods 49,50 falls in between those from Sato et al. 48 and McSkimin and Andreatch 51 . Modifying the higher-order elastic properties of the diamond is another advancement over ref. 27 .  Fig. 13 Schematic of evolution of the yield surface f ðp; s; ε p ; ε t p ðτÞÞ ¼ 0 until it reaches the fixed surface of perfect plasticity φðp; sÞ ¼ 0 in a "five-dimensional" space of deviatoric stresses s at fixed pressure p. The initial yield surface and the fixed surface of perfect plasticity φðp; sÞ ¼ 0 are isotropic and pictured here as circles with their center at O. Two other yield surfaces depend on plastic strain ε p at the current time t and the entire plastic strain history ε t p ðτÞ before time t. These surfaces acquire strain-induced anisotropy, which is depicted by the shifted centers of the surfaces O 1 and O 2 . However, when the yield surface reaches the fixed surface of perfect plasticity φðp; sÞ ¼ 0, which is isotropic and plastic strain-and plastic strain history-independent, moves along it at further loading, the material deforms as perfectly plastic, isotropic with the fixed surface of perfect plasticity φðp; sÞ ¼ 0 V.I. Levitas et al.