# Computational prediction of ideal strength for a material

Zixun Wang1, Xingyu Wang2, Xianqi Song1, Xinxin Zhang3,
Hanyu Liu1,*, Miao Zhang2,*

1 _Key Laboratory of Material Simulation Methods and Software of Ministry of
Education & State Key Laboratory of Superhard Materials, College of Physics,
Jilin University, Changchun 130012, China_

_2 Department of Physics, School of Sciences, Beihua University, Jilin 132013,
China_

3 _Department of Science, Shenyang University of Chemical Technology, Shenyang
110142, China_

*Corresponding author:

_hanyuliu@jlu.edu.cn_(Hanyu Liu);

_zhangmiaolmc@126.com_(Miao Zhang)

## Abstract

The ideal strength is crucial for predicting material behavior under extreme
conditions, which can provide insights into material limits, guide design and
engineer for enhanced performance and durability. In this work, we present a
method within an allows for the estimation of tensile, shear, and indentation
strengths in any crystallographic direction or plane. We have examined the
strain-stress relationships of several well-known structures and compared our
findings with previous work, demonstrating the effectiveness of our approach.
Moreover, we performed extensive investigations into the indentation strength
of hexagonal WC, $\beta$-SiC, and MgAl2O4. The current study uncovers the
modes of structural deformation and the underlying atomistic mechanisms. The
insights gained from this study have significant implications for the further
exploration and design of superhard materials.

Keywords: First-principles calculations; Ideal strength; Tensile strength;
Shear strength; Indentation strength.

## Introduction

The ideal elasticity strength is the minimum stress to plastically deform an
infinite dislocation-free crystal [1], which gives an upper bound on the
strength of a perfect crystal. The ideal strength could be classified as
tensile, shear, and indentation strength according to the type of external
force posted on the material. It is an important parameter to evaluate the
shear and indentation strengths are believed to be closely related to the
Vickers hardness of an ideal crystal. Therefore, it is of technical and
scientific interest to investigate the ideal elasticity strength of a crystal.
In industrial practice, stress-strain relationships are commonly employed to
determine ideal strength [2], where various strains are applied to the crystal
and the corresponding stress is measured. Theoretical simulation of the ideal
strength under complex loading conditions by using stress-strain relationship
can also provide fundamental insights into the evolution of chemical bonding
and physical property changes of a materials, which is helpful for elucidating
and predicting anomalous behaviors of material at extreme conditions [2, 3, 4,
5].

Great efforts had been dedicated to the accurate calculation of the ideal
strength since the mid of 1980s[1]. Initially, some experimentally proposed
semi-empirical models were used. It is commonly believed that the tensile and
shear strength of materials are directly proportional to their Vickers
hardness. However, the results are always inaccurate. Recent advances in
computational physics made it possible to ab-initially calculate the ideal
strength, and the ideal elasticity strength of tremendous materials has been
precisely calculated within the framework of electronic structure[1, 6, 7, 8,
9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29]. Indentation is a method used in experiments to measure the hardness and
strength of materials. A more complex calculation of indentation strength has
been performed in the framework of DFT [16, 27, 28]. Energy and force can be
calculated faster and more accurately using the DFT approach. Therefore, this
calculation can be applied in more scenarios. The calculations of the elastic
strength have been extensively applied to an increasing number of materials,
solving a wide range of engineering and physical problems [1, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29].
Li _et al_. investigated boron-rich tungsten borides, and demonstrated novel
ultrahard properties with stress-strain behaviors revealing diverse and
anomalous patterns. The profound influence of boron concentration on bonding
configurations and deformation modes in WBn (n = 2, 3, 4) results in distinct
stress responses and unexpected variations in indentation strength,
highlighting the unique deformation mechanisms and the need to explore
unconventional structure-property relations in ultrahard materials.[21]Lu et
al. designed as new-generation superhard materials, transition-metal light-
element compounds have been hindered by indentation strain softening, limiting
their intrinsic hardness below 40 GPa. Moreover, the calculations revealed
that hP4-WN and hP6-WN2 tungsten nitrides exhibit extraordinary strain
stiffening, leading to enhanced indentation strengths exceeding 40 GPa,
opening possibilities for nontraditional superhard materials.[24]

In this work, we provided a comprehensive exposition of the underlying
principles of our newly proposed methodology on calculating the tensile, shear
and indentation strengths based on Vienna Ab initio Simulation Package (VASP)
code, accompanied by comparative analysis with other approaches to establish
its robustness. We selected hexagonal WC, $\beta$-SiC, and MgAl2O4 for an in-
depth investigation of ideal strength. Additionally, we examined the
indentation strength of hexagonal WC, $\beta$-SiC, and MgAl2O4.This
comprehensive exploration has afforded us a deeper comprehension of their
mechanical characteristics.

## Computational Methods

Mathematically, the stress tensors are the second-order symmetric tensor with
nine components (three positive stress componentsand six shear stress
components), in which only six of them are independent. They also can be
represented by a 3×3 matrix, where each matrix element corresponds to a stress
component. To calculate the stress-strain relationship, it should strain the
crystal in a series of incremental simple tensile or shear, and then relax the
lattice as a function of the strain to obtain the total energy _E_ tot.
According to DFT, _E_ tot can be obtained by

${E_{tot}}=\sum\limits_{i}{\left.{\left\langle{{\psi_{i}}\left|{\left.{-\frac{1}{2}{\nabla^{2}}}\right|{\psi_{i}}}\right.}\right.}\right\rangle}+\int{{V_{ext}}}(r){d^{3}}r+{E_{H}}(\rho)+{E_{xc}}(\rho)+{U_{II}}$

where _$\Psi$ i_ represents wavefunction, _V ext_ is external potential, _E
H_(_$\rho$_) is the Hartree energy, _E xc_(_$\rho$_) is exchange and
correlation energy, and _U II_ is the ion-ion interaction. _$\rho$(r)_ is the
density of the electron gas, which can be calculated by

$\rho(r){\text{ = 2}}{\sum\limits_{i}{\left|{{\psi_{i}}}\right|}^{2}}$

Based on Hellmann-Feynman force theorem, each stress component is the
derivative of the energy with respect to the strain

${\tau_{mn}}=-\frac{1}{V}\frac{{\partial{E_{tot}}}}{{\partial{\varepsilon_{mn}}}}=-\left.{\left\langle{{\psi_{i}}\left|{\left.{\frac{{\partial\widehat{H}}}{{\partial{\varepsilon_{mn}}}}}\right|{\psi_{i}}}\right.}\right.}\right\rangle=-\frac{{\partial{U_{II}}}}{{\partial{\varepsilon_{mn}}}}-\int{\frac{{\partial{V_{ext}}(r)}}{{\partial{\varepsilon_{mn}}}}\rho({\varepsilon_{mn}}){d^{3}}{\varepsilon_{mn}}}$

where H represents the Hamiltonian, and _$\varepsilon$ mn_ represents strain.

According to the above method, VASP can simply determine the six independent
stress tensors _$\tau$ xx_, _$\tau$ yy_, _$\tau$ zz_, _$\tau$ xy_, _$\tau$
xz_, _and $\tau$yz_ by a structural optimizations, if a strain is induced in a
well optimized crystal lattice by modifying the lattice parameters, which are
represented by a lattice matrix in VASP code. For instance, if calculating the
tensile stress along [100] direction in fcc lattice, one just needs to add the
tensile strain along the [100] direction (or x-axis of crystal lattice) by
changing the corresponding lattice matrix element, and then optimize the
structure by fixing the _F_ xx component of the optimization matrix _F_ mn (m,
n = x, y, z) in VASP code without optimization. Then we can obtain the tensile
stress of _$\tau$ xx_ which corresponds to the studied tensile stress
_$\tau$_[100]. Considering that VASP can straightforwardly compute the stress
tensor of _$\tau$ xx_, _$\tau$ yy_, _$\tau$ zz_, if the tensile strain was
added to the coordinate axis of crystal lattice. Therefore, we rotate the
studied tensile direction [_uvw_] to one of coordinate axis (here, we choose
_x_ -axis). In this way, after structural optimization by setting _F_ xx to
zero, we can obtained the tensile stress of _$\tau$ xx _which corresponds to
the studied tensile stress _$\tau$_[_uvw_]. By gradually applying a series of
continuous strain along the tensile direction and ultimately deriving the
corresponding tensile stress, we can obtain the tensile stress-strain
relationship.

It is obvious that during the tensile strength calculations, one important
step is to rotate the studied tensile direction to the _x_ -axis. Here, take
the orthorhombic lattice as an example to introduce how to rotate the studied
tensile direction to _x_ -axis. Using Cartesian coordinate system, the actual
stretching direction is OA as shown in Fig. 1(a). It can be achieved by
rotating -$\alpha$ degree (here counterclockwise rotation is set to be
positive direction) around the _z_ -axis (see Fig. 1b), then rotate $\beta$
degree around the _y_ -axis.

Fig. 1 Rotation of tensile direction in an orthohormbic lattice. (a) tensile
direction before rotation. The red vector denotes the studied tensile
direction, the green line is its projection on _xy_ -plane, $\alpha$ is the
angle of projection to the x-axis, $\beta$ is the angle of tensile direction
to xy-plane. (b) The lattice after tensile direction rotating -$\alpha$ around
_z_ -axis, in which the tensile direction is located in _xz_ -plane.

Shear strength is more complicated than tensile strength. Shear stress
calculations need to rotate the studied shear direction [_uvw_] to z-axis and
reorient the normal vector of the shear plane to align with _x_ -axis.
Consequently, by maintaining Fxz as zero during structural optimization, we
can calculate the shear stress _$\tau$ xz,_ which is corresponding to the
studied shear stress _$\tau$_[_uvw_]. Here, to give a more vivid illustration
on the rotation process of shear plane and shear direction, orthorhombic
lattice was also taken as an example. As it is can be seen in Fig. 2(a), the
sheared plane is characterized by its normal vector OA and the shear direction
OB lies within this plane. Firstly, the normal vector OA was rotated -$\alpha$
degree around _z_ -axis [see Fig. 2(b)], and then $\beta$ degree around _y_
-axis [see Fig. 2(c)]. Subsequently, rotate OA $\gamma$ degree around the
x-axis[Fig. 2(d)]. Finally, the normal vector is align with the x-axis, and
the shear stress direction aligns with the z-axis.

Fig. 2 Rotation of shear plane and direction in an orthohormbic lattice. (a)
The orthorhombic lattice with shear plane and direction marked out. The yellow
vector represents shear direction, and the red vector represents normal vector
of the shear plane. The green line is the projection on xy-plane. $\alpha$ is
the angle of projection to the x-axis, and $\beta$ is the angle of tensile
direction to xy-plane. (b) The lattice after rotating -$\alpha$ around z-axis
to make normal direction in xz-plane. (c) The lattice after rotating $\beta$
around y-axis to make tensile direction in x-axis. (d) The lattice after
rotating $\gamma$ around y-axis to make shear direction in z-axis.

The indentation strength is a more complex case, where in actual compression
the material receives not only tensile stress perpendicular to the plane, but
also in-plane shear stress. Experiments have shown that there is a close
relationship between nano-scale indentation and the geometry of the indenter
[16, 22, 23, 29]. When using the first-principles method to calculate
compression, we assume that the material is uniformly deformed within a small
volume. And the deformation is related to the shape of the indenter. Such a
uniformly deformed lattice is sufficient in the study of the breakage of
interactions on the atomic microscopic scale. In the first-principles
calculations, we can simplify the situation to a uniform deformed infinite and
perfect lattice under a uniform stress field[30]. We assume that a sharp
indenter squeezes the material by the way shown in the Fig. 3(a). The tensile
stress $\tau$_xx_ and shear stress $\tau$_xz_ received by the material under a
sharp indenter always maintain the following relationship: $\tau$_xz_
=$\tau$_xx_ tan$\theta$, where $\theta$ is determined by the shape of the
indenter, as shown in Fig. 3(b).

Fig. 3 Microscopic indentation. (a) A sharp indentation acts on a material,
and 2$\theta$ is the angle of indenter. (b) The microscopic lattice beneath
the indenter.

Based on this relationship, we only need to calculate the shear stress by
using the method of calculating shear stress proposed in this work. Note that
since we consider the crystal is an infinitely large perfect lattice with
uniform deformation, we must keep the angle of the deformation constant in
both directions. This can be achieved by modifying the optimization program in
VASP code. For a blunt indenter, as long as we know that the direction of the
tangent at each point of contact surface between the indenter and the
material, we can use the approach of a sharp indenter to calculate the
indentation stress at each point of the blunt indenter, as shown in Fig.
4[30].

Fig. 4 Blunt indenter.

## Results and discussion

To thoroughly validate the accuracy and reliability of our approach, we
conducted an extensive analysis involving diamond, Si, hexagonal AlN, and
$\alpha$-Al2O3. The tensile, shear and indentation strengths were calculated
separately by setting the relative strain at each step to 0.01 times the
elements of the lattice matrix. We used the VASP package with the Perdew-
Burke-Ernzerhof (PBE) version of the generalized gradient approximation (GGA)
for the exchange-correlation functional and a plane-wave basis set [31, 32,
33, 34]. For the integration in the Brillouin zone, we set the k-points to 6 ×
6 × 6 for all the calculated material. The cutoff energy is set to 600 eV,
which is sufficient to make the total energy converge to an accuracy of less
than 0.001 eV/atom. The paramount objective of our investigation is to
establish the veracity of our proposed methodology. In pursuit of this goal,
we performed a meticulous comparison of our computed results against those
obtained using established methods. This validation process was meticulously
conducted under conditions consistent with prior research, ensuring a fair and
direct comparison.

The crystallographic configurations of the chosen materials are presented in
Table S1, while the outcomes of the tensile, shear, and indentation are
elaborated in Tables 1, respectively. The significance of these tables lies in
their role as a tangible demonstration of the effectiveness of our method
across various loading scenarios. In the case of diamond, our calculations
yielded the weakest tensile and shear directions consistent with previous
research. The weakest tensile direction is [111], and the weakest shear
direction is [111] <$11\overline{2}$>. Additionally, our Vickers hardness
value is 95.7 GPa, which closely aligns with the prior result of 97.6 GPa,
showcasing minimal deviation. This conforms to the characteristics of diamond
as a superhard material. For cases Si, the weakest tensile and shear
directions are [111] and [111] <$11\overline{2}$>, respectively. And its
indentation strength is 6.1 GPa. Additionally, the weakest tensile direction
for hexagonal AlN is [$\overline{1}2\overline{1}0$], and the weakest shear
direction is [$\overline{1}2\overline{1}0$] <$10\overline{1}0$>. The
indentation strength is 11.7 GPa. Moreover, the weakest tensile and shear
directions are [$10\overline{1}0$] and [$10\overline{1}0$] <0001> for
$\alpha$-Al2O3. The indentation strength is 28.9 GPa. Our computed results
also closely match those from the previous work. The specific stress-strain
curves for tension, shear, and indentation of these four materials are
presented in Figures S1-4, respectively.

Table 1: Ideal tensile, shear and indentation strength in three directions for
diamond, Si, hex-AlN, and _$\alpha$_ -Al2O3, where hex denotes hexagonal
phase.

Structure | Tensile strength (GPa) | Shear strength (GPa) | Indentation strength (GPa)
---|---|---|---
| Direction | This work | Other | Direction | This work | Other | Direction | This work | Other
diamond | [100] | 225.7 |  | [100] <001> | 140.4 |  | [111] <11$\overline{2}$> | 95.7 | 97.6e
| [110] | 126.4 |  | [110] <001> | 98.5 |  |  |  |
| [111] | 92.4 | 95a | [111] <11$\overline{2}$> | 93.7 | 93.0a |  |  |
Si | [100] | 86.8 |  | [100] <001> | 12.1 |  | [111] <11$\overline{2}$> | 6.1 |
| [110] | 29.5 |  | [001] <100> | 15.0 |  |  |  |
| [111] | 21.1 | 22b | [110] <001> | 7.5 | 6.8b |  |  |
Hex-AlN | [0001] | 38.1 | 39c | [10$\overline{1}$0] <0001> | 20.2 | 20.0c | [10$\overline{1}$0] <$\overline{1}$2$\overline{1}$0> | 11.7 |
| [10$\overline{1}$0] | 38.0 | 38c | [$\overline{1}$2$\overline{1}$0] <0001> | 20.0 | 22.5c |  |  |
| [$\overline{1}$2$\overline{1}$0] | 32.1 | 34c | [10$\overline{1}$0] <$\overline{1}$2$\overline{1}$0> | 19.6 | 19.9c |  |  |
| [0001] | 57.3 | 58.3d | [10$\overline{1}$0] <0001> | 16.4 |  |  |  |
$\alpha$-Al2O3 | [10$\overline{1}$0] | 26.3 | 26.3d | [$\overline{1}$2$\overline{1}$0] <0001> | 14.5 | 14.4d | [$\overline{1}$2$\overline{1}$0] <0001> | 28.9 |
| [$\overline{1}$2$\overline{1}$0] | 36.5 | 36.5d | [10$\overline{1}$0] <$\overline{1}$2$\overline{1}$0> | 18.7 |  |  |  |

  * •

a Ref [12] b Ref [18] c Ref [10] d Ref [11] e Ref [16]

We subsequently placed particular emphasis on an in-depth investigation of the
ideal strengths of three distinct materials: WC, $\beta$-SiC,and MgAl2O4. We
would like to calculate the Vickers hardness, so the indenter we consider to
be a Vickers indenter with an angle of 136°between the opposing surfaces.

Hexagonal tungsten carbide (WC) is crucial for cutting-edge industries like
manufacturing, aerospace, and electronics, thanks to its exceptional
properties. Its unique combination of hardness, strength, and thermal
stability significantly enhances overall performance and efficiency. We
extensively undertook an investigation into the ideal tensile strength, ideal
shear strength, and indentation strength of hexagonal tungsten carbide (WC).
According to academic standards, Fig.5(a) illustrates several distinctive
crystallographic orientations of the WC crystal, which predominantly guided
our calculations in this study. We conducted tensile calculations to explore
these properties along the [0001], [$10\overline{1}0$] and
[$\overline{1}2\overline{1}0]$ crystallographic directions. The tensile
strengths were found to be 101.3 GPa, 72.6 GPa, and 63.5 GPa in the [0001],
[$10\overline{1}0$] and [$\overline{1}2\overline{1}0]$ orientation,
respectively. Fig.5(b) visually presents the stress-strain curves for the
investigated directions, highlighting [$\overline{1}2\overline{1}0]$ as the
least resistant tensile orientation within WC.

We next undertook an in-depth examination involving shear calculations of
properties along the [$\overline{1}2\overline{1}0]$<0001>,
[$10\overline{1}0$]<0001>, and [$\overline{1}2\overline{1}0$]
<$10\overline{1}0$> crystallographic orientations. The shear strength along
the [$\overline{1}2\overline{1}0]$<0001> direction emerged at an impressive
41.6 GPa, while the [$10\overline{1}0$]<0001> orientation showcased a
noteworthy tensile strength of 33.9 GPa. Additionally, the
[$\overline{1}2\overline{1}0$] <$10\overline{1}0$> orientation unveiled a
distinctive tensile strength of 18.7 GPa. These significant findings are shown
in Fig.5(c), elegantly capturing the essence of stress-strain profiles for
each respective orientation. Remarkably, [$\overline{1}2\overline{1}0$]
<$10\overline{1}0$> represents the weakest shear direction of WC. Therefore,
we selected the weakest shear direction as the orientation for its
indentation.

A comprehensive analysis of the indentation process was conducted and
visualized in Fig.5(d). Intriguingly, as the strain surpassed 0.27 times the
crystal length in the aforementioned direction, the chemical bonds between
tungsten and carbon atoms suffered rupture due to increased atomic separation.
At this point, the indentation strength of WC is 36.8GPa. Its stress-strain
curve is displayed in Fig.5 (e). This phenomenon is attributed to the
exceeding distance between the atoms, leading to bond fracture. Our study
provides valuable insights into the mechanical response of hexagonal WC under
various loading conditions and crystallographic orientations, shedding light
on the intricate interplay between atomic bonds and mechanical deformation.
The calculation results indicate that WC is a relatively hard material with
the potential to become a superhard material.

Fig. 5 (a)Structure and the crystallographic directions of WC (b) The tensile
stress-strain curve. (c)The shear stress-strain curve. (d) Structural changes
during the indentation process and fracture of chemical bonds. (e)The
indentation stress-strain curve.

At equilibrium condition, _$\beta$_ -SiC adopts the cubic structure in _F_ -43
_m_ symmetry [Fig. 6(a)], where each Si atom is connected to four C atoms,
forming a strong _sp_ 3 covalent bond unit with Si-C bond lengths of 1.89 Å.
We first set out to determine the easy cleavage plane of _$\beta$_-SiC by
simulating its ideal tensile strength along several high-symmetry directions.
The calculated results shown in Fig. 6(b) clearly indicate that the weakest
peak tensile stress occurs in the [111] direction in _$\beta$_ -SiC, giving
rise the (111) easy cleavage planes. The result is 45 GPa in <111> direction
for _$\beta$_ -SiC and the Si-C bonds soften considerably under the tensile
strain long before the bond breaking. Since a material’s hardness is primarily
determined by its ability to resist shear deformation, we examine shear stress
responses of _$\beta$_-SiC in its major low-index planes to establish key
benchmarks. We show in Fig.6(c) the stress responses under pure and
indentation shear strains along the (111) [11$\overset{\overline{}}{2}$]
direction. The two sets of shear stresses reach peak values 29.7 and 19.7 GPa,
respectively. At the indentation strain ($\varepsilon$ = 0.36), the stress has
only decreased from the peak value. The peak pure shear stress occurs at the
strain $\varepsilon$ = 0.22 but the interlayer _sp_ 3 bonding persists up to
$\varepsilon$ = 0.32 where significant buckling still exists and the Si-C bond
length increased to 2.076 Å. The corresponding structural snapshot shows
significant buckling in the (111) planes, indicating a strong _sp_ 3 bonding
character. The structural failure occurs at $\varepsilon$ = 0.38, where the
interlayer bonds in the (111) [11$\overset{\overline{}}{2}$] direction break
up and the layers become essentially flat, signaling a dominant _sp_ 2 bonding
character. Here we notice that all the Si-C bonds responsible for the
structural failure also break simultaneously under shear strains.The
calculation results show that _$\beta$_ -SiC possesses favorable mechanical
properties.

FIG.6. (a) The crystal structure of _$\beta$_ -SiC at equilibrium (top) and
projected onto the (111) plane (bottom). (b) Calculated stress responses under
tensile strains along various high-symmetry directions of _$\beta$_-SiC. (c)
Calculated stress responses under pure shear and indentation shear strains in
the (111) plane along the [11$\overset{\overline{}}{2}$] shear direction. (d)
Bond-length changes under the (111) [11$\overset{\overline{}}{2}$] pure shear
and indentation strains. Corresponding snapshots of selected bonding
configurations at key deformation points before and after the large drop of
stress on each stress-strain curve under pure shear and indentation strains
are presented in (e) and (f), respectively.

MgAl2O4 holds significant importance across diverse fields due to its
remarkable combination of properties. As a transparent ceramic, it finds
applications in optical systems, lasers, and transparent armor. Its high-
temperature stability, electrical insulating properties, and chemical
resistance make it valuable in aerospace, electronics, and corrosion-resistant
coatings. Moreover, spinel’s biocompatibility has fueled its adoption in
medical implants. Its versatile attributes and extensive utility underscore
spinel’s role as a multifunctional material with widespread impact.

Our investigation centered around the characterization of spinel MgAl2O4 in
terms of its ideal tensile strength, ideal shear strength, and indentation
strength. To adhere to academic conventions, we adhered to Fig.7(a), which
visually illustrates distinct crystallographic orientations within the MgAl2O4
crystal. These orientations served as the foundation for our computations in
this study. Tensile calculations were carried out along the [100], [110], and
[111] crystallographic directions, yielding corresponding strengths of 89.2
GPa, 28.28 GPa, and 28.3 GPa, respectively. The stress-strain curves for these
orientations are depicted in Fig.7(b), with the [111] orientation revealing
the least resistance to tensile forces within MgAl2O4.

Our subsequent focus was on an exhaustive exploration of shear calculations
along the $[100]$<001>, $[110$]<001>, and [111] <$11\overline{2}$>
crystallographic orientations. Notably, the shear strength along
the$[100]$<001> direction registered an impressive 29.7 GPa, while the
$[110$]<001> orientation displayed a remarkable tensile strength of 31.0 GPa.
Furthermore, the [111] <$11\overline{2}$> orientation exhibited a distinctive
tensile strength of 14.2 GPa. Fig6 (c) elegantly encapsulates these pivotal
findings, showcasing stress-strain profiles for each orientation.
Interestingly, the [111] <$11\overline{2}$> orientation emerged as the weakest
shear direction within MgAl2O4, thereby being chosen for indentation analysis.

Our comprehensive analysis of the indentation process is vividly presented in
Fig. 7(d). Intriguingly, when the strain surpassed 0.46 times the crystal
length in the aforementioned direction, the increased atomic separation led to
the rupture of chemical bonds. At this juncture, the indentation strength of
MgAl2O4 was determined to be 22.2 GPa, as depicted in the stress-strain curve
presented in Fig.7(e). This phenomenon was attributed to the exceeding atomic
distance, which induced bond fracture. Our study thus yields valuable insights
into the mechanical behavior of hexagonal MgAl2O4 under varied loading
conditions and crystallographic orientations, unraveling the intricate
interplay between atomic bonds and mechanical deformation. Its hardness
enables it to be applied in more diverse applications.

Fig. 7 (a)Structure and the crystallographic directions of MgAl2O4 (b) The
tensile stress-strain curve. (c)The shear stress-strain curve. (d) Structural
changes during the indentation process and fracture of chemical bonds. (e)The
indentation stress-strain curve.

## Conclusions

In summary, we introduced a novel methodology for computing the ideal strength
of materials based on first-principles calculations. The approach enabled the
determination of tensile, shear, and indentation strengths in arbitrary
crystallographic directions, facilitating a comprehensive understanding of
material response under different loading conditions. Through a meticulous
comparison with prior works, we demonstrated the accuracy and reliability of
our proposed method. Furthermore, our investigation delved into the specific
cases of hexagonal WC, cubic SiC, and spinel MgAl2O4, shedding light on their
ideal strengths and the structural transformations that occur under applied
stress. Based on our computational results, the Vickers hardness of WC reaches
36.8 GPa, which is closed to 40 GPa. Therefore, it could not be classified as
a superhard material. While WC has the potential to become a superhard
material through various methods, such as nanocrystallization, high-pressure
synthesis, and alloying. Meanwhile, cubic SiC and MgAl2O4 exhibit Vickers
hardness within the range of 20-30 GPa, rendering them suitable for diverse
applications in various fields. This detailed analysis not only enriched our
understanding of these materials’ mechanical properties but also highlighted
the versatility of our methodology in exploring a diverse range of crystalline
structures. The findings presented in this study have broad implications for
material design and engineering, offering a rigorous foundation for predicting
and optimizing the mechanical behaviors of various materials. As we continue
to refine and expand our computational techniques, the insights gained from
this research will contribute to advancements across multiple disciplines,
driving innovation in materials science and engineering.

## Acknowledgments

This work is supported by the Natural Science Foundation of China (Grant
No.12074138), and the Jilin Provincial Science and Technology Development
Joint Fund Project (Grant No. YDZJ202201ZYTS581).

## References

  * [1] D. Roundy, C. R. Krenn, Marvin L. Cohen, and J. W. Morris.  Ideal shear strengths of fcc aluminum and copper.  Phys. Rev. Lett., 82:2713–2716, Mar 1999.
  * [2] Weizhong Han, Ling Huang, Shigenobu Ogata, Hajime Kimizuka, Zhao-Chun Yang, Christopher R. Weinberger, Qing-Jie Li, Boyu Liu, Xi xiang Zhang, Ju Li, Evan Ma, and Zhi wei Shan.  From “smaller is stronger” to “size-independent strength plateau”: Towards measuring the ideal strength of iron.  Advanced Materials, 27, 2015.
  * [3] D.F. Bahr, D.E. Kramer, and W.W. Gerberich.  Non-linear deformation mechanisms during nanoindentation.  Acta Materialia, 46(10):3605–3617, 1998.
  * [4] A. Kelly and N.H. Macmillan.  Strong Solids.  Monographs on the physics and chemistry of materials. Clarendon Press, 1986.
  * [5] Qian Yu, Zhi-Wei Shan, Ju Li, Xiaoxu Huang, Lin Xiao, Jun Sun, and Evan Ma.  Strong crystal size effect on deformation twinning.  Nature, 463(7279):335–338, Jan 2010.
  * [6] M. de Jong, W. Chen, T. Angsten, A. Jain, R. Notestine, A. Gamst, M. Sluiter, C.K. Ande, S. van der Zwaag, J.J. Plata, C. Toher, S. Curtarolo, G. Ceder, K.A. Persson, and M. Asta.  Charting the complete elastic properties of inorganic crystalline compounds.  Scientific Data, 2, March 2015.
  * [7] Vadim V. Brazhkin, Alexander G. Lyapin, and Russell J. Hemley.  Harder than diamond: Dreams and reality.  Philosophical Magazine A, 82(2):231–253, 2002.
  * [8] Y. Sun and E. Kaxiras.  Slip energy barriers in aluminium and implications for ductile versus brittle behaviour.  Phil. Mag. A, 75:1117, 1997.
  * [9] Yi Zhang, Hong Sun, and Changfeng Chen.  Superhard cubic ${\mathrm{b}\mathrm{c}}_{2}\mathrm{N}$ compared to diamond.  Phys. Rev. Lett., 93:195504, Nov 2004.
  * [10] R. F. Zhang, S. H. Sheng, and S. Veprek.  Mechanism of the $b3$ to $b1$ transformation in cubic aln under uniaxial stress.  Phys. Rev. B, 76:075208, Aug 2007.
  * [11] Z. Zhang, Z. H. Fu, R. F. Zhang, D. Legut, and H. B. Guo.  Anomalous mechanical strengths and shear deformation paths of al2o3 polymorphs with high ionicity.  RSC Adv., 6:12885–12892, 2016.
  * [12] H. Chacham and Leonard Kleinman.  Instabilities in diamond under high shear stress.  Phys. Rev. Lett., 85:4904–4907, Dec 2000.
  * [13] Seung-Hoon Jhi, Steven G. Louie, Marvin L. Cohen, and J. W. Morris.  Mechanical instability and ideal shear strength of transition metal carbides and nitrides.  Phys. Rev. Lett., 87:075503, Jul 2001.
  * [14] D. M. Clatterbuck, C. R. Krenn, Marvin L. Cohen, and J. W. Morris.  Phonon instabilities and the ideal strength of aluminum.  Phys. Rev. Lett., 91:135501, Sep 2003.
  * [15] X. Blase, Philippe Gillet, A. San Miguel, and P. Mélinon.  Exceptional ideal strength of carbon clathrates.  Phys. Rev. Lett., 92:215505, May 2004.
  * [16] Zicheng Pan, Hong Sun, and Changfeng Chen.  Indenter-angle-sensitive fracture modes and stress response at incipient plasticity.  Phys. Rev. B, 79:104102, Mar 2009.
  * [17] Zicheng Pan, Hong Sun, and Changfeng Chen.  Indenter-angle-sensitive fracture modes and stress response at incipient plasticity.  Phys. Rev. B, 79:104102, Mar 2009.
  * [18] David Roundy and Marvin L. Cohen.  Ideal strength of diamond, si, and ge.  Phys. Rev. B, 64:212103, Nov 2001.
  * [19] Yinwei Li, Jian Hao, Hanyu Liu, Siyu Lu, and John S. Tse.  High-energy density and superhard nitrogen-rich b-n compounds.  Phys. Rev. Lett., 115:105502, Sep 2015.
  * [20] Bingjia Xu, Yingxiao Mu, Zhu Mao, Zongliang Xie, Haozhong Wu, Yi Zhang, Chongjun Jin, Zhenguo Chi, Siwei Liu, Jiarui Xu, Yuan-Chun Wu, Po-Yen Lu, Alan Lien, and Martin R. Bryce.  Achieving remarkable mechanochromism and white-light emission with thermally activated delayed fluorescence through the molecular heredity principle.  Chem. Sci., 7:2201–2206, 2016.
  * [21] Quan Li, Dan Zhou, Weitao Zheng, Yanming Ma, and Changfeng Chen.  Anomalous stress response of ultrahard ${\mathrm{wb}}_{n}$ compounds.  Phys. Rev. Lett., 115:185502, Oct 2015.
  * [22] Shigenobu Ogata, Ju Li, and Sidney Yip.  Ideal pure shear strength of aluminum and copper.  Science, 298(5594):807–811, 2002.
  * [23] Shigenobu Ogata, Ju Li, Naoto Hirosaki, Yoji Shibutani, and Sidney Yip.  Ideal shear strain of metals and ceramics.  Phys. Rev. B, 70:104104, Sep 2004.
  * [24] Cheng Lu, Quan Li, Yanming Ma, and Changfeng Chen.  Extraordinary indentation strain stiffening produces superhard tungsten nitrides.  Phys. Rev. Lett., 119:115503, Sep 2017.
  * [25] Zhijun Shi, Wei Shao, Lixiang Rao, Xiaolei Xing, Yefei Zhou, Xuejun Ren, and Qingxiang Yang.  Pressure-induced structure, elasticity, intrinsic hardness and ideal strength of tetragonal c4n.  Phys. Chem. Chem. Phys., 24:5171–5184, 2022.
  * [26] Lei Chen, Junlian Xu, Meiguang Zhang, Gangtai Zhang, Yaru Zhao, Yun Zhang, and Congcong Ren.  Theoretical study on novel orthorhombic ternary monocarbides m0.5re0.5c (m=v, nb, ta) from first-principles calculations.  Ceramics International, 46(15):24624–24634, 2020.
  * [27] Lin Wen and Aiyuan Hu.  Effect of biquadratic exchange and anisotropy on the critical temperature of antiferromagnet.  Acta Physica Sinica, 69:107501, 2020.
  * [28] Cheng Lu and Changfeng Chen.  Indentation-strain stiffening in tungsten nitrides: Mechanisms and implications.  Phys. Rev. Mater., 4:043402, Apr 2020.
  * [29] A. Gouldstone, H.-J. Koh, K.-Y. Zeng, A.E. Giannakopoulos, and S. Suresh.  Discrete and continuous deformation during nanoindentation of thin films.  Acta Materialia, 48(9):2277–2295, 2000.
  * [30] H. Bei, E. P. George, J. L. Hay, and G. M. Pharr.  Influence of indenter tip geometry on elastic deformation during nanoindentation.  Phys. Rev. Lett., 95:045501, Jul 2005.
  * [31] J Ihm, A Zunger, and M L Cohen.  Momentum-space formalism for the total energy of solids.  Journal of Physics C: Solid State Physics, 12(21):4409, nov 1979\.
  * [32] D. M. Ceperley and B. J. Alder.  Ground state of the electron gas by a stochastic method.  Phys. Rev. Lett., 45:566–569, Aug 1980.
  * [33] John P. Perdew, Kieron Burke, and Matthias Ernzerhof.  Generalized gradient approximation made simple.  Phys. Rev. Lett., 77:3865–3868, Oct 1996.
  * [34] G. Kresse and J. Furthmüller.  Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.  Phys. Rev. B, 54:11169–11186, Oct 1996.


# STABILITY DIAGRAM OF COLLIDING BEAMS

X. Buffat  EPFL  Lausanne; CERN  Geneva  Switzerland
W. Herr  N. Mounet  T. Pieloni  CERN  Geneva  Switzerland

###### Abstract

The effect of the beam–beam interactions on the stability of impedance mode is
discussed. The detuning is evaluated by the means of single particle tracking
in arbitrarily complex collision configurations, including lattice non-
linearities, and used to numerically evaluate the dispersion integral. This
approach also allows the effect of non-Gaussian distributions to be
considered. Distributions modified by the action of external noise are
discussed.

##  1 Introduction

The stability of impedance driven modes is usually ensured by non-linearities
via Landau damping. The strength of the Landau damping is represented in a
stability diagram, defining the area of complex tune shifts which can be kept
stable [1]. The computation of the stability diagram is a crucial element for
the evaluation of the non-linearity required for the beam stability. It can be
done analytically in many simple cases, in particular considering lattice non-
linearities and a limited number of beam–beam interactions. Facing the
complexity of the LHC collision scheme, a numerical approach seems more
appropriate. In particular, the dispersion integral is solved numerically
using the amplitude detuning provided by single particle tracking, thus
allowing the computation of stability diagrams in arbitrarily complex cases.

Following the same approach, the effect of non-Gaussian distribution functions
can be introduced in the computation. Previous studies have shown that the
distribution function plays a crucial role in the stability diagram [2] and,
being usually poorly known, greatly diminishes the predicting power of such a
consideration. The numerical evaluation of the stability diagram allows us to
go further in the understanding of this effect by using a non-analytical
distribution function. In particular, coloured external noise is known to
enhance diffusion of resonant particles [3], greatly modifying the
distribution and therefore the stability diagram.

##  2 Stability of multibunch coherent beam–beam modes

The method to derive the stability of impedance driven multibunch modes by
placing their tune shifts in the stability diagram derived using the lattice
non-linearities, as described in [4], considers bunches with identical
detuning. In the LHC, this assumption is no longer valid once beam–beam
effects are considered. Indeed, the asymmetric layout of the interaction
points as well as the asymmetric filling scheme lead to a variety of bunches
having significantly different detuning. Theoretical treatment of the beam
stability in such configurations is currently lacking. It is however possible
to assess these cases using multiparticle tracking simulation [5], at the
expense of large computational power.

The beam–beam interactions do not only introduce bunch dependent detuning,
they also change the nature of the coherent modes. The Landau damping of
beam–beam modes in the single bunch regime is addressed in [6], the extension
to multibunch coherent beam–beam mode is, however, not trivial. Preliminary
studies by the means of multiparticle tracking simulation are presented in
[5]. Such an approach is well suited to studying the LHC in its full
complexity, however it is very demanding in terms of computing power. It is
therefore interesting to consider simplified cases. In this paper, we address
the stability of single bunch impedance modes with detuning from the lattice
and beam–beam interactions, in other words, multibunch effects from the
impedance as well as beam–beam coherent mode are neglected. These drastic
assumptions are motivated by the observation, in the LHC, of single bunch
instabilities, while operating in the multibunch regime. Figure 1 shows such
an instability at the end of the squeeze during operation of the LHC in 2012.
In particular, the measurement of the beam oscillation amplitude provided by
the BBQ system indicates a coherent instability on Beam 1. The observation of
bunch intensities indicates that only one bunch lost its intensity in an
abnormal way, with respect to the other bunches, suggesting that only this
bunch had undergone the instability.

(a) BBQ (b) FBCT

Figure 1: Observation of an instability during luminosity production in the
LHC (fill 2644). The beam oscillation amplitude shows a coherent excitation in
the horizontal plane of Beam 1. The bunch by bunch intensity of both beams
shows that all bunches (in blue) except one (in red) are stable.

##  3 Numerical evaluation of stability diagrams

Considering a case without coupling, the stability diagram of each plane is
obtained by solving the dispersion relation for a given detuning
$q(J_{x},J_{y})$ and distribution function $\psi(J_{x},J_{y})$ where $J_{x}$
and $J_{y}$ are the unperturbed actions in each plane,

$\frac{-1}{\Delta
Q_{i}}=\iint_{0}^{\infty}\frac{J_{i}\frac{\mathrm{d}\psi}{\mathrm{d}J_{i}}}{Q-q(J_{x},J_{y})}\mathrm{d}J_{x}\mathrm{d}J_{y},~{}Q\in\mathbb{R},~{}i=1,2.$
(1)

The $\Delta Q_{i}$ found for different values of $Q$ are the tune shifts at
the limit of stability, therefore, they define an area in which the tune
shifts are stable. In cases where the denominator is strictly non-zero, the
tune shift is purely real, which indicates the absence of Landau damping. In
the relevant cases, the integral possesses poles, and can be evaluated using
various techniques. Whereas it is difficult to obtain analytically in complex
configurations, the detuning is easily obtainable numerically. The dispersion
integral can then be evaluated by standard numerical techniques, in our case
by adding a vanishing complex part to the denominator.

###  3.1 Detuning

Figure 2: Illustration of a typical tune footprint with octupoles (positive
polarity) and long-range beam–beam interactions at the end of the squeeze in
the LHC. The tracking result is often distorted by the presence of resonances;
using a fine mesh and a repairing algorithm allows us to get a footprint
suitable for the computation of stability diagrams.

The detuning is obtained numerically from tracking simulation with MAD-X [7,
8, 9]. Particles with different actions are tracked for 1024 turns from which
the tunes can be evaluated, the result is represented in the tune diagram as a
tune footprint (Fig. 2). A significant number of particles is required to
reduce the noise in the evaluation of the integrals. In particular, the effect
of resonant driving term on the tracking can have a significant effect on the
evaluation of the footprint. This effect is reduced by an automatic removal of
faulty points and interpolation from the closest non-faulty points. As
illustrated in Fig. 2, this technique does not allow us to fully remove
distortions due to strong non-linearities, in particular from coupling
resonances. In the cases of practical interests, however, the effect on the
footprint, and consequently on the stability diagram, is tolerable. The
integrity of the footprint is nevertheless systematically checked before
deriving a stability diagram.

###  3.2 Benchmarking

Figure 3 shows a good agreement between stability diagrams computed
analytically and numerically in the case of octupolar detuning using a
Gaussian distribution.

Practical applications of this code for the LHC can be found in [10].

Figure 3: Comparison of analytical and numerical derivation of a stability
diagram with LHC octupoles powered with -100 A for a 4 TeV beam with a
normalized emittance of 2 $\mu$m.

##  4 The effect of external noise

External coloured noise enhances the diffusion of resonant particles. In the
presence of amplitude detuning, this results in a depletion of certain parts
of the distribution. This effect is illustrated with a simple model, using a
multiparticle code that tracks particles through a lattice with amplitude
detuning described by

$\displaystyle Q_{x}$ $\displaystyle=$ $\displaystyle Q_{x,0}+a\cdot
J_{x}+b\cdot J_{y},$ (2) $\displaystyle Q_{y}$ $\displaystyle=$ $\displaystyle
Q_{y,0}+b\cdot J_{x}+a\cdot J_{y}.$ (3)

The noise is modelled with a sinusoidal excitation with finite correlation
time. The relative difference to the initial distribution in action space
after $2\cdot 10^{4}$ turns is shown in Fig. 4LABEL:sub@fig-distr_noise. The
depletion of the distribution at the position of resonant particles is
visible. This effect was already studied for different purposes in [11]. The
measurable effect on the distribution in real space is very small (Fig.
4LABEL:sub@fig-distr_real), the effect on the stability diagram is, however,
significant, as illustrated by the comparison of the stability diagrams
obtained using the initial and perturbed distributions (Fig. 4LABEL:sub@fig-
sdiag_noise).

(a) Relative difference to initial distribution in action space after $2\cdot
10^{4}$ turns

(b) (Un)perturbed stability diagrams

(c) Real space bunch profile

Figure 4: Stability diagram derived from a distribution perturbed by external
noise in the presence of amplitude detuning $a=1.64\cdot 10^{-4},b=-1.16\cdot
10^{-4}$. The noise has an amplitude of $10^{-4}\sigma_{x^{\prime}}$ and a
frequency of $Q_{x,0}+0.05$ with a correlation time of 100 turns.

The effect illustrated by this simple model can be simulated in more realistic
scenarios. In particular, the COMBI code [12], extended with an impedance
model, was used to model the action of external noise on a beam stabilized by
both amplitude detuning and a transverse feedback. We consider a study case,
where a single bunch is rendered unstable by negative chromaticity. This study
case intends to mimic a situation similar to the LHC at the end of the squeeze
where the beams have to be stabilized by high damper gain and high octupole
current, while relaxing significantly the need for computational power due to
the large number of bunches in this case. The results are therefore not meant
to be compared to measurements, but rather to illustrate an effect that should
be studied in more realistic cases in the future. The upper dot in Fig.
5LABEL:sub@fig-study_case_sdiag shows the coherent tune shift of the most
unstable mode, with a large positive imaginary part. The imaginary part is
brought down to the lower dot using a transverse feedback. The line represents
the stability diagram provided by an octupolar amplitude detuning. In this
configuration, this analysis suggests that both the damper and the octupoles
are required to stabilize the beam, which is confirmed by tracking (Fig.
5LABEL:sub@fig-study_case_rise), showing the oscillation amplitude as a
function of time for the different scenarios. The stable configuration can,
however, be rendered unstable by a wideband noise, i.e. a random kick is
applied every turn, the kick being constant over the bunch length. Indeed, the
response to the wideband noise is more important at the mode frequencies as it
cannot be entirely damped by a transverse feedback with finite gain. Wideband
noise in this configuration has a similar effect on the diffusion as coloured
noise, and therefore can compromise the stability from Landau damping, as
described above. In particular, in Fig. 6, the simulation of the stable case
for different noise amplitude shows that the latency time before the
instability takes off depends on the noise amplitude.

(a) Stability diagram and tune shifts

(b) Tracking without external noise

Figure 5: Analysis of a single bunch with intensity $1.5\cdot 10^{11}$ in the
LHC, chromaticity of $-10$ units, stabilized with a transverse damper with
gain $\sim 700$ turns and amplitude detuning $a=-2.05\cdot 10^{-4},b=1.45\cdot
10^{-4}$.

(a) Tracking in the stable configuration with different noise amplitude (b)
Latency time with a noise amplitude of $5\cdot 10^{-5}~{}\sigma_{x^{\prime}}$

Figure 6: Tracking simulations in the stable configuration (Fig. 5) rendered
unstable by wideband noise. The linear fit is done on the first part, and the
exponential fit during the instability, which indicates a latency of $2.8\cdot
10^{5}$ turns before an instability of rise time $2.7\cdot 10^{4}$ turns.

The extension of this study case to more realistic configurations including
multibunch impedance and beam–beam interactions promises great challenges from
the computational point of view, but is however not out of reach with current
resources. In particular, one could expect that other sources of diffusion may
counterbalance the effect described. Nevertheless, this analysis suggests that
the external noise could play a critical role in configurations where both a
transverse feedback and amplitude detuning are required to stabilize the
beams.

##  5 Conclusion

Single bunch stability of beams colliding in arbitrary complex configurations
can be investigated by numerically solving the dispersion integral, using
amplitude detuning derived from single particle tracking simulations. A new
code, based on MAD-X output, has been developed to assess the stability of the
LHC beams in any operational configuration. This approach is however not
suited to studying the stability of multibunch coherent modes.

The effect of distributions perturbed by external noise on the stability
diagram has been investigated using the code mentioned above, and revealed
very strong effects for modifications of the distributions well below the
sensitivity of current measurements. It has been shown using multiparticle
tracking codes that, indeed, a beam well stabilized by both amplitude detuning
and transverse damper could be rendered unstable by the introduction of
wideband noise.

##  6 Acknowledgments

The authors would like to acknowledge S. White and E. Métral for fruitful
discussions.

## References

  * [1] A.W. Chao, Physics of Collective Beams Instabilities in High Energy Accelerators, (John Wiley and Sons, New York, 1993).
  * [2] E. Métral $et~{}al$., “Stability Diagram for Landau Damping with a Beam Collimated at an Arbitrary Number of Sigmas,” CERN-AB-2004-019-ABP (2004).
  * [3] A. Bazzani $et~{}al$., “Diffusion in Hamiltonian Systems Driven by Harmonic Noise,” J. Phys. A: Math. Gen. 31 (1998).
  * [4] J.P. Koutchouk and F. Ruggiero, “A Summary on Landau Octupoles for the LHC,” LHC Project Note 163 (1998).
  * [5] S. White $et~{}al$., “Beam–Beam and Impedance,” these proceedings.
  * [6] Y.Alexahin, “On the Landau Damping and Decoherence of Transverse Dipole Oscillations in Colliding Beams,” CERN SL-96-064 AP (1996).
  * [7] http://mad.web.cern.ch/mad/.
  * [8] W. Herr, “Particle Tracking With MAD-X Including LHC Beam–Beam Interactions,” LHC Project Note 344 (2004).
  * [9] W. Herr, “Tune Shift and Spread Due to Short and Long Range Beam–Beam Interactions in the LHC,” CERN SL/90-06 (AP), LHC Note No. 119 (1990).
  * [10] X. Buffat $et~{}al$., “Operational Considerations on the Stability of Colliding Beams,” these proceedings.
  * [11] A. Bazzani $et~{}al$., “Effect of Colored Noise on the Betatronic Motion: A Possible Mechanism for Slow Extraction,” Nonlinear and Stochastic Beam Dynamics in Accelerators - A Challenge to Theoretical and Computational Physics, Lüneburg (1997).
  * [12] T. Pieloni, “A Study of Beam–Beam Effects in Hadron Colliders with a Large Number of Bunches,” EPFL Thèse No 4211 (2008).