  * Keywords: few-layer graphene, rhombohedral stacking, Bernal stacking, deformation, stability

# Relative Stability of Bernal and Rhombohedral Stackings in Trilayer Graphene
under Distortions

R Guerrero-Avilés1 , M Pelc2,3, F R Geisenhof3, T Weitz4, and A Ayuela1 1
Material Physics Center CFM-MPC, Donostia International Physics Center (DIPC),
Paseo Manuel Lardizabal 4-5, 20018 Donostia-San Sebastián, Spain  2 Institute
of Physics, Nicolaus Copernicus University in Toruń, Grudziadzka 5, 87-100
Toruń, Poland  3 Physics of Nanosystems, Department of Physics, Ludwig-
Maximilians-Universität München, Amalienstrasse 54, 80799 Munich, Germany  4
I. Physikalisches Institut, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany
<EMAIL_ADDRESS>
###### Abstract

Stackings in graphene have a pivotal role in properties to be discussed in the
future, as seen in the recently found superconductivity of twisted bilayer
graphene. Beyond bilayer graphene, the stacking order of multilayer graphene
can be rhombohedral, which shows flat bands near the Fermi level that are
associated with interesting phenomena, such as tunable conducting surface
states expected to exhibit spontaneous quantum Hall effect, surface
superconductivity, and even topological order.However, the difficulty in
exploring rhombohedral graphenes is that in experiments, the alternative,
hexagonal stacking is the most commonly found geometry and has been considered
the most stable configuration for many years.Here we reexamine this stability
issue in line with current ongoing studies in various laboratories. We
conducted a detailed investigation of the relative stability of trilayer
graphene stackings and showed how delicate this subject is. These few-layer
graphenes appear to have two basic stackings with similar energies. The
rhombohedral and Bernal stackings are selected using not only compressions but
anisotropic in-plane distortions. Furthermore, switching between stable
stackings is more clearly induced by deformations such as shear and breaking
of the symmetries between graphene sublattices, which can be accessed during
selective synthesis approaches. We seek a guide on how to better control – by
preserving and changing - the stackings in multilayer graphene samples.

## 1\. Introduction

Graphene remains a key two-dimensional material to still successfully reveal
intriguing properties when stacked into just a few layers. As an example,
bilayer graphene in twisted form has been found to exhibit superconductivity
properties [1]. Large graphene samples are obtained experimentally by
mechanical exfoliation of graphite, so they are supposed to favor the Bernal
stacking, which for graphite has been considered more stable than the
rhombohedral stacking [2]. However, this stability is far from being fully
settled. Firstly, recent calculations showed that graphite is energetically
more stable in the rhombohedral stacking, so it is suggested that crystal
growing conditions under high temperatures and pressures in a geological scale
can stand behind graphite in the Bernal stacking [3]. Secondly, experimental
and synthesizing techniques in 2D materials have been so improved over the
years that graphene is being investigated nowadays with a defined number of
layers [4, 5]. Research works that focused on the same exfoliated sample of
trilayer graphene (TLG) simultaneously found regions with the Bernal and
rhombohedral stacking [6, 7, 8, 9, 10]. Rhombohedral graphenes are today being
synthesized with a revival in a research race to show unique properties mostly
due to their flat bands [11]. The stability between the Bernal and
rhombohedral stacking in few layer graphenes therefore stands as a crucial
question to be discussed.

Techniques are currently being developed on how to obtain, preserve or change
stacking in trilayer graphene (TLG). The phase change between rhombohedral
(ABC) and Bernal (ABA) stacking can be induced by external driving forces. The
ABC stacking was experimentally found when graphene layers were grown on SiC
substrates, and it was stable up to a temperature of around 1,200 C∘ when it
reverts to the ABA stacking [12, 13]. The transition from ABA to ABC stacking
can be caused by a scanning tunneling microscope (STM) tip on highly ordered
pyrolytic graphite (HOPG) samples due to the small barriers involved ($\sim
1.0$ mev/atom), as recently reported [7, 14]. The transition between ABC and
ABA stacking was also achieved using triazine decoration, which causes a large
energy difference to induce the stacking transition [15]. The stacking in few-
layer graphenes was stabilized mechanically by strain, using electric fields
and even doping [16, 17, 18]. In order to fabricate high-quality encapsulated
ABC trilayer - hBN devices used for quantum transport investigations, it has
been shown that often the ABC trilayer converts to ABA during van-der-Waals
transfer, even though ABA regions had been carefully removed prior to stamping
[19]. The rhombohedral to Bernal transition has been theoretically studied by
modeling the continuous sliding deformation of a side layer and by
parameterizing results for two layers, calculations that were performed within
density functional theory (DFT) using local density approximation (LDA) [20,
21]. Although the van-der-Waals energies per atom are small, the total bonding
energies are sizeable when the sample areas are considered, total energies
that are large but needed to keep a certain phase locked in a metastable
stacking after being synthesized. These findings suggest that not only the
relative stability between these TLG stackings should be studied, but that
their dependence on the lattice deformations deserves further study.

In this work, we address the stability between rhombohedral and Bernal
stacking in trilayer graphene by calculations performed within density
functional theory (DFT) using explicitly van der Waals functionals. We study
different deformations that induced the energetic exchange between these
stackings by looking at the relative stability. We first show that as the TLG
ground state, the rhombohedral stacking is slightly more stable than the
Bernal one. Then we perform several deformations in the quest to induce the
transition from rhombohedral stacking to the Bernal one. In our simulations we
find that the Bernal stacking is energetically more stable than the
rhombohedral one under some distortions based on the planar expansion of the
lattice parameters and their perpendicular compression. Furthermore, we show
that under distortions such as shear and the atomic shifts in graphene
sublattices, the structural anisotropy helps to stabilize the ABA stackings in
certain directions. These mechanisms point out that the stacking transition is
achieved by breaking the sublattices symmetries of layer-layer interactions
involved in each stacking. These findings indicate that different practical
realizations are possible, such as depositing samples in sandwiches,
substrates, and molecule decoration, where stacking transitions may take
place.

## 2\. TLG stability and electronic structures

Figure 1(a) shows the relaxed unit cells of trilayer graphene with the Bernal
stacking on the left and the rhombohedral stacking on the right. Each stacking
has a relaxed intralayer atomic distance of $a_{c-c}=1.43$Å, and the
interlayer distance of $d=3.55$Å in agreement with experiments [22]. The
relative stability between Bernal and rhombohedral stackings is next analyzed
by looking at the difference between their total energies. After geometry
relaxations, the energy differences per atom are calculated as $\Delta
E=(E_{B}-E_{r})/N$, where $E_{B}$ and $E_{r}$ are the total energies of Bernal
and rhombohedral stacking respectively, and $N=6$ is being the number of atoms
in the unit cells. For pristine TLG our former results show a total energy
difference of 0.079 meV/atom between both stackings. We have performed a
detailed analysis of converging parameters to the accuracy required here, as
shown in supplementary data. It seems theoretically established that the
Bernal stacking is more stable than rhombohedral stacking in the literature,
but this statement is based on inconsistent results [3]. In fact, we find that
in TLG the rhombohedral stacking is slightly more favorable than the Bernal
one.

Figure 1: (a) (Left) Side views of the units cells and (right) band structures
for the Bernal and rhombohedral stacking of trilayer graphene (TLG). (b)
Energy differences between Bernal and rhombohedral TLG versus in-plane
homogeneous strains $\delta$ with compression ($\delta<$0) and stretching
($\delta>$0) in lattice parameters. (c) Energy differences for uniaxial
deformations along zigzag and armchair directions. Blue-cyan and red-orange
colors refer to the strain values when the Bernal and rhombohedral stackings
are more stable, respectively. Energy differences are also given with strains
assuming a constant area per atom, shown in light red. The insets in panels
(b) and (c) indicate the homogeneous and uniaxial deformations in-plane,
respectively. The TLG rhombohedral stacking is more stable without surface
strain, and a transition to Bernal stacking is shown under expansions with
$\delta\geq 1\%$.

The energy differences although small are in agreement with enthalpy
experiments between hexagonal and rhombohedral graphene [18, 23]. Kinetic
conditions such as temperatures and pressures seem to be the key in obtaining
the Bernal stacking. All atoms have to shift the entire layer in order to pass
from ABA to ABC stacking, and the other way around. We propose that these
processes have to be considered in order to compare with the room temperature
energy taking into account the total contact area. This means that there is a
minimum area in the flakes as to obtain a stacking stable at room temperature.
Depending on the value of stacking energy times area, there is a blocking
temperature above which the separate layers can slide. In other words there is
also a minimum area at a given temperature so that the flake stacking remain
stable to be observed in experiments. An implication is that stackings can be
different depending on the synthesis methods using exfoliation and CVD growth.
We compare the electronic band structures of the Bernal and rhombohedral
stackings along the $\Gamma$KM path near the Fermi energy, as shown in Figure
1(a). At the K point the Bernal stacking bands are displaying superposed a
massive and massless electronic behavior, and the rhombohedral electronic
structure shows massive electrons only; these two types of dispersion bands
for each stacking are in agreement with Ref. [24]. The Bernal stacking can be
seen as a superposition of single layer graphene bands and Bernal bilayer
graphene bands; however, the graphene-like bands are slightly shifted up
indicating the source of some intrinsic electric field between the layers to
be added to external fields [25]. For the rhombohedral stacking there is a
parabolic behavior with a flat region around the K point, source of many
interesting properties, such as tunable conducting surface states [11]
expected to exhibit spontaneous quantum Hall effect [26], surface
superconductivity [27], and even topological order [28]. The massless band in
the Bernal stacking now has acquired mass in the rhombohedral geometry. The
conduction and valence bands are degenerated and slightly displaced from the K
point in the KM direction. At lower energies, the rhombohedral stacking
presents split bands what justifies its stability found in the calculations. A
simple image can allow us to explain this stability. The Bernal stacking has
directly connected three atoms between layers so that there are states at zero
energy. However, the rhombohedral case has two pairs of atoms between layers
so the corresponding electronic states are split from the zero energy, a fact
that could be somehow linked to its higher stability.

## 3\. In-plane deformations

We now study different in-plane deformations of the lattice vectors. Each
deformation is characterized by a stretch factor $\lambda$ that operates in
the lattice vectors on different ways. First, we consider a homogeneous (or
isotropic) deformation of lattice vectors described by
$\vec{a}_{1,2}=\lambda(a_{x}\hat{i}+a_{y}\hat{j})$. In this sense, having
$\lambda$=1 (or strain $\delta=0\%$) shows the previous relaxed lattice
vectors, while a value of $\lambda=1.05$ indicates an expansion strain of
5$\%$ ($\lambda=0.95$ is a compression of $5\%$). This causes a proportional
deformation of the lattice vectors keeping the interlayer distance constant.
Figure 1(b) shows the total energy differences per atom between Bernal and
rhombohedral stacking for the in-plane deformations explained above. Positive
values of the energy differences (colored in red) indicate that the
rhombohedral stacking is more stable than the Bernal one. On the contrary,
negative values for the energy differences (colored in blue) refer to the case
of the Bernal stacking being more stable. We show that although under
homogeneous deformation the Bernal stacking becomes more favorable on
expansion, the energy difference decreases for positive strains till it
changes sign for $\sim 1\%$. At this point the Bernal stacking becomes more
favourable.

Next, we consider anisotropic deformations with uniaxial strain without and
with constant area, described by the lattice vectors $\vec{a}_{1,2}=\lambda
a_{x}\hat{i}+a_{y}\hat{j}$ ($\vec{a}_{1,2}=a_{x}\hat{x}+\lambda a_{y}\hat{y}$)
and $\vec{a}_{1,2}=\lambda a_{x}\hat{i}+a_{y}/(\lambda)\hat{j}$ along
$\hat{x}$ ($\hat{y}$) direction, respectively. The energy difference induced
by the anisotropic deformation considering a constant area is slightly
decreasing and remains almost constant in the whole range of considered
positive or negative strains in the order of 5 %, as shown in Figure 1 (c).
Furthermore, this figure includes the energy differences with the anisotropic
uniaxial deformations, and we find that the Bernal stacking becomes more
stable during stretching. The values of energy in the strain range $\pm 5\%$
show smaller response compared with those of homogeneous deformation, 2 versus
5 tenths of meV. The rhombohedral stacking with anisotropic strains is more
stable for strain values till $\delta\sim 2.2\%$, to be compared with the
Bernal stacking with isotropic strains that becomes more stable under small
strains $\delta\sim 1.1\%$ .

Figure 2 summarizes the trends found between the Bernal and rhombohedral
energy differences versus the in-plane deformations in TLG. There is an
asymmetry in strains that induces that the most stable rhombohedral stacking
switches into the Bernal stacking. The strain values when this transition is
taking place are within the experimental range of reversible strains attained
for samples in substrates, already found for a single layer of graphene [29].
In fact, strain values of beyond 1% could be artificially applied to graphene
heterostructures [30] and strain up to 0.3% has been seen in graphene
encapsulated in hBN [31]. Thus detailed care must then be taken in experiments
when depositing exfoliated graphenes on substrates, and when using patterned
contacts [10, 32], so that strain is induced by top contacts while edge
contacts can help in reducing strain [33].

Figure 2:  Summary of energy differences with respect to in-plane deformations
$\delta_{x}$ and $\delta_{y}$. Note the anisotropy of the rhombohedral-Bernal
stability versus in-plane deformations.

## 4\. Out of plane deformations

We further study the relative stability order induced by the strain $\epsilon$
in the interlayer distance by stretching and compression between layers, as
shown in Figure. 3. Although the energies of both stackings are almost
degenerated at zero strain, $\Delta E$ under stretching shows that the
rhombohedral continues being more stable even for an interlayer distance up to
$d=4.26$Å ($\epsilon\leq 20\%$). For a compression at $d=3.11$Å
($\epsilon=-2.5\%$), the $\Delta E$ values show that the Bernal stacking
becomes more favorable, in agreement with another current calculation using
empirically fitted vdW Grimme functional [34]. Small compressions induced by
processes during transfer and protective layers will be enough to end into the
ABA stacking. While we are not aware of direct measurements of compressions
during h-BN encapsulation or subsequent cleaning by squeezing out residues
[35], one can anticipate that compression of a trilayer can take place during
these processing steps.

Figure 3:  Out of plane deformations: total energy differences versus out-of-
plane strain $\epsilon$ related to the interlayer distance.

The next question is to determine how under a perpendicular pressure also the
atomic in-plane layers are modified. It is noteworthy that under constant
volume a perpendicular pressure implies that the layers are expanding on the
plane. For instance assuming a constant volume, a value of $\epsilon=-3.5\%$
corresponds to a homogenous expansion of $\delta=1.88\%$. The critical
$\delta$ value reinforces that the transition to Bernal stacking comes coupled
to the in-plane deformations, as shown in Figure 1(b) and (c) above. These
results indicate that the TLG stackings have to be analyzed in the sample once
they are covered and isolated in the operating device whenever is possible to
be certain that the rhombohedral stacking has not changed into Bernal.

## 5\. Shear strain deformations

To assess the role of symmetry breaking we analyze the shear strain
deformations. We consider the shear deformations along the perpendicular layer
direction using the azimuthal $\phi$ and polar $\theta$ angles, as shown in
the inset of Figure 4 and the movie included in supplementary data. The
question is what $\theta$ values must be chosen before entering into different
stackings from the rhombohedral and Bernal ones, i.e. how close the structures
with shear can be seen as deformations of the original ABA and ABC stackings.
The shear structures remain close to the pristine ABC and ABA stackings for
the values with $\theta<15^{\circ}$, as shown in supplementary data, where
shears beyond this critical angle show layers that can hardly be seen anymore
as small deformations of the Bernal and rhombohedral stackings.

Figure 4 shows energy differences between Bernal and rhombohedral stackings
for the shears with $\theta=10^{\circ},5^{\circ}$ changing $\phi$ angles. The
shear for $\theta=10^{\circ}$ also provokes that the interlayer distance _d_
decreases from the value in pristine lattices so that $\epsilon=-1.5\%$;
however, this $\epsilon$ value is still well below the critical distance to
get the ABA stacking with perpendicular pressure. So that when obtaining
Bernal stackings is not just correlated with the distance but better with the
shear angle. Note that for $\theta\neq 0^{\circ}$ the shear is displacing the
layers in the in-plane with respect to the perfect stacking as cones, shown in
the upper panel, so that the overlap of the carbon pz orbitals between layers
change related not only to distances but also to angles. In fact the
$\epsilon$ value for $\theta=5^{\circ}$ ($\epsilon\sim 0.3\%$) just decreases
slightly the interlayer distance, being larger than the critical distance to
obtain ABA stacking under pressure.

Figure 4: Stability energy difference between Bernal and rhombohedral
stackings versus shear deformation, which is indicated by angles $\theta$ and
$\phi$ shown in the upper panel. The inset shows the shear $\phi$ angles for
which the ABA and ABC stackings are more stable.

On the one hand, the energies of Bernal stacking exhibit a $60^{\circ}$
periodicity in $\phi$, while increasing $\theta$ makes the energy profile just
wider. That means that these deformations destabilize Bernal stacking
minimally up to $\phi=30^{\circ}$. On the other hand, the energies for
rhombohedral stacking does not show the same periodicity along $\phi$ angle,
because their periodicity (120∘ in $\phi$) is twice than one found for Bernal
stacking. Thus we find a region where the rhombohedral stacking is more
destabilized than the Bernal one.

In Figure 4 we show the energy difference per atom between the two stackings
under shear deformations. The curve for $\theta=5^{\circ}$ shows the $\Delta
E$ in the range of $\sim\pm 0.1meV/atom$ decreasing as $\phi$ values increase.
The case of $\theta=10^{\circ}$ where the $\Delta E$ are in the range of
$\sim\pm 6$ tenths of meV/atom much above the hundreths of meV/atom found for
previous studied deformations; in fact they become an order of magnitude
larger. For values of $\phi$ less than $30^{\circ}$ the rhombohedral stacking
still remains as more favorable, and for the values of $\phi>30^{\circ}$, we
find that the Bernal stacking is more stable than the rhombohedral one. Taking
into account the symmetry, we find that the Bernal stacking becomes favorable
for three fold directions of $\phi=60^{\circ}$ as shown in the inset included
in Figure 4. Interestingly, the ABA stackings for $\theta=5^{\circ}$ are
observed even for very small interlayer distance, values with very small
$\epsilon$. Even for a test case with $\epsilon=0$, i.e with no change in
interlayer distance, the energies remain in the order of tenths of meV. In
general, considering that those $\Delta$E values are larger than previous
considered deformation cases, the shear deformation of the stackings promote
the Bernal stacking in certain directions, so that the area of stacking
regions can even be reduced by a factor ten to be stable using shear
deformations. These findings are in line with previous experimental results in
which shear is induced by contacts [10], by other encapsulating layered
graphenes within BN 111Along zigzag and armchair direction the ABC and ABA
stackings were favored, respectively [16], and by explicitly putting shear to
the samples [36].

## 6\. Perpendicular sublattice distortions

Up to now all the investigated deformations affected both graphene sublattices
within the same layer in the same way. We have then shown that shear
introduced mostly anisotropy in plane. For the sake of completeness we now
asses the role of deformations out of plane by looking at the anisotropy
induced by the sublattices sites. We study the displacement up and down for
each of the sublattices in both stackings. We characterize how much the same
perpendicular shift affects the relative stability of each stacking. Here, we
apply a perpendicular displacement to A or B sublattice type nodes, shifts
that are breaking the sublattice layer symmetry. We shift an atom of a
particular sublattice in the direction perpendicular to the layer plane by
less than the layer compression to get ABA stacking, and compare how much the
total energy of each stacking is changed.

Figure 5:  (a) Models for the Bernal and rhombohedral stacking including the
number of interlayer next nearest neighbors of the middle layer atoms, to be
discussed in text. (b) Comparison of relative energies with respect the ground
state of the perpendicular displacements of the A${}_{2}^{\uparrow}$ and
B${}_{2}^{\uparrow}$ graphene sublattices in the middle layer. The subscript
denotes layer number, and the superscript arrow indicates the corresponding up
and down shift of the sublattice atom.

First we note that in each stacking, atoms of both side layers have the same
number of interlayer neighbors. The energy response for atom shifts in side
layers is very small, it does not exceed 0.01 meV. Thus we focus on the middle
layer analyzing the coordination numbers of its atoms. Figure 5 (a) shows the
Bernal and rhombohedral stacking schemes indicating the nearest neighbours in
next layers. The vertical lines are linking the middle layer atoms to the
first neighbors belonging to the side layers and shaded prims with numbers
mark the second interlayer neighbors. In both stackings the A and B type nodes
have in total 2 first neighbors and 18 second neighbors belonging to the side
layers. However, the number of neighbors for each sublattice calculated
separately is different. We can expect different responses when breaking the
symmetry by shifting A or B atoms. In the case of Bernal stacking the numbers
of interlayer neighbors for a sublattice site are the same with top and bottom
layer. Thus the up or down shifts will be symmetric. In the rhombohedral
stacking the numbers of neighbors for each sublattice site are “antisymmetric”
i.e. shifting the A atom up is equivalent to shift the B type atom down and
vice versa. For the rhombohedral geometry the number of cases to be discussed
is reduced to two: (i) $A_{2}^{\uparrow}$ equivalent to $B_{2}^{\downarrow}$
and (ii) $A_{2}^{\downarrow}$ equivalent to $B_{2}^{\uparrow}$, where
subscripts refer to the layer number, and arrows indicate the perpendicular
displacement.

Figure 5(b) shows the comparison of energy response difference defined as
$E_{shift}-E_{0}$, where for each stacking in the two types of displacements
$E_{0}$ is the ground state energy, and $E_{shift}$ is the total energy when
the atoms are shifted. We shift the atoms by 1% of interlayer distance, i.e.
0.0335Å. In all the considered cases the total energy in each stacking
increases by $\approx$ 0.5 meV. The energy response due to each sublattice
shift remains nearly the same for the ABA stacking, and for the ABC stacking
it changes more depending on the shift direction. For ABC stacking, moving the
A atom up - equivalent to shift the B atom down - causes different energy
response in the two directions. In the case of rhombohedral stacking, the
energy increase is higher by 0.07 meV. As a result the Bernal stacking becomes
more stable. The situation is just opposite when moving the A atom down
(equivalent to shift the B atom up) - in this case the Bernal stacking is more
destabilized, and the rhombohedral geometry remains favourable. The
perpendicular sublattice shift studied here is a lattice deformation that
causes an energy difference between stackings about 0.07 meV. This finding
points out that breaking the sublattice symmetry is another way to stabilize
the Bernal and rhombohedral stackings.

These energy differences due to breaking of sublattice symmetry can be
explained by having the TLG stackings interpreted again in terms of trimers
and dimers between layers. The total energy can be written in two parts as a
sum of eigenvalues and a sum of interatomic potentials with other atoms.
Looking at the total contribution of interatomic potentials the second
neighbours atoms in the nearest layers are held responsible for these energy
differences. This energy contribution adds to the difference in the sum of
eigenvalues, already commented, in which the levels for the case of Bernal
stacking split from the zero energy when the atoms are shifted up or down.
These calculations raise further intringuing experiments regarding how the
sublattice symmetry of stackings could be considerably broken by substrates
and adatoms, or even better by having molecules and nanoparticles conforming
these graphene heterostructures [37, 38].

## 7\. Conclusions

The aim of the present paper was to determine when the graphene stacking
changes from rhombohedral to Bernal due to small lattice deformations. Under
compression the rhombohedral stacking suffers a phase transformation to Bernal
with strains $\epsilon<2.5\%$. For even smaller compressions, shear
deformations and shifts breaking graphene sublattice symmetries induce an
anisotropy that stabilizes the Bernal stacking. These findings provide
insights into the role of substrate-associated strains when graphene layers
are integrated into devices, so that the stacking order and consequently their
ultimate electronic properties are modified. This work would be of interest in
relevant technological areas such as patterning contacts, and encapsulating
graphene flakes between other materials. Our results are then asking for
further experiments looking at the role of directional shears along and the
adsorption of adatoms and molecules on graphene stackings.

This research was funded by the Spanish Ministry of Science and Innovation
(grants no. PID2019-105488GB-I00), the Gobierno Vasco UPV/EHU (project no.
IT-1246-19), and the European Commission NRG-STORAGE project (project no. GA
870114). This research was conducted in the scope of the Transnational Common
Laboratory (LTC) “Aquitaine-Euskadi Network in Green Concrete and Cement-based
Materials.”

## Methods

We performed our study within density functional theory (DFT) using the
ab-$initio$ simulation package VASP [39, 40, 41]. This method based on plane
waves is applied using a well converged kinetic energy cut-off of 700 Ry. In
the energy range that we are interested, the mesh grid on the reciprocal space
is key to get consistent results (see supplementary data), a
288$\times$288$\times$1 centered at the $\Gamma$ point is used. For the
dispersive interactions between layers, we perform our calculations using the
Van der Waals functional vdW-DF2 [42]. The electronic convergence was
performed with $10^{-7}$ eV. We relax the cell shape and keep the volume
constant, and we fix the ions in the $xy$ plane while allowing them to move in
the $z$ direction. Tests performed with other functionals, as well as detailed
convergence test in k-points, and electronic and ionic relaxations are
included in supplementary data. Note that a fine k-mesh is required to be
applied with the number of sampling k points in the x and y directions being a
multiple of 3 to explicitly include the points K, K´and a region around them,
regions in the reciprocal space where the low energy physics of multilayer
graphene occurs.

## Supplementary data

The following files are available free of charge.

  1. $\bullet$

Convergence tests to properly simulate few-layer graphene. Tests using Several
vdW Functionals. Farthest shear deformations in Bernal and rhombohedral
stackings.

  2. $\bullet$

Movie to visually follow the shear deformations of trilayer graphene
stackings.

This material is available free of charge via the Internet at
https://iopscience.iop.org/.

## References

## References

  * [1] Cao Y, Fatemi V, Fang S, Watanabe K, Taniguchi T, Kaxiras E and Jarillo-Herrero P 2018 Nature 556 43–50
  * [2] Tománek D and Louie S G 1988 Physical Review B 37(14) 8327–36
  * [3] Xiao R, Tasnadi F, Koepernik K, Venderbos J, Richter M and Taut M 2011 Physical Review B 84 165404–18
  * [4] Sugawara K, Yamamura N, Matsuda K, Norimatsu W, Kusunoki M, Sato T and Takahashi T 2018 NPG Asia Materials 10 e466–e466
  * [5] Park J, An H, Choi D C, Hussain S, Song W, An K S, Lee W J, Lee N, Lee W G and Jung J 2016 Nanoscale 8(30) 14633–42
  * [6] Lui C H, Li Z, Chen Z, Klimov P V, Brus L E and Heinz T F 2011 Nano Letters 11 164–9
  * [7] Xu P, Yang Y, Qi D, Barber S D, Ackerman M L, Schoelz J K, Bothwell T B, Barraza-Lopez S, Bellaiche L and Thibado P M 2012 Applied Physics Letters 100 201601–4
  * [8] Cong C, Yu T, Sato K, Shang J, Saito R, Dresselhaus G F and Dresselhaus M S 2011 ACS Nano 5 8760–8
  * [9] Warner J H, Mukai M and Kirkland A I 2012 ACS Nano 6 5680–6
  * [10] Geisenhof F R, Winterer F, Wakolbinger S, Gokus T D, Durmaz Y C, Priesack D, Lenz J, Keilmann F, Watanabe K, Taniguchi T et al. 2019 ACS Applied Nano Materials 2 6067–75
  * [11] Shi Y, Xu S, Yang Y, Slizovskiy S, Morozov S V, Son S K, Ozdemir S, Mullan C, Barrier J, Yin J et al. 2020 Nature 584 210–4
  * [12] Ferrari A C, Meyer J C, Scardaci V, Casiraghi C, Lazzeri M, Mauri F, Piscanec S, Jiang D, Novoselov K S, Roth S and Geim A K 2006 Physical Review Letters 97(18) 187401–6
  * [13] Graf D, Molitor F, Ensslin K, Stampfer C, Jungen A, Hierold C and Wirtz L 2007 Nano Letters 7 238–42
  * [14] Yankowitz M, Joel I, Wang J, Birdwell A G, Chen Y A, Watanabe K, Taniguchi T, Jacquod P, San-Jose P, Jarillo-Herrero P et al. 2014 Nature Materials 13 786–9
  * [15] Zhang W, Yan J, Chen C H, Lei L, Kuo J L, Shen Z and Li L J 2013 Nature Communications 4 2074–81
  * [16] Yang Y, Zou Y C, Woods C R, Shi Y, Yin J, Xu S, Ozdemir S, Taniguchi T, Watanabe K, Geim A K et al. 2019 Nano Letters 19 8526–32
  * [17] Gao Z, Wang S, Berry J, Zhang Q, Gebhardt J, Parkin W M, Avila J, Yi H, Chen C, Hurtado-Parra S et al. 2020 Nature Communications 11 1–10
  * [18] Li H, Utama M I B, Wang S, Zhao W, Zhao S, Xiao X, Jiang Y, Jiang L, Taniguchi T, Watanabe K et al. 2020 Nano Letters 20 3106–12
  * [19] Chen G, Sharpe A L, Gallagher P, Rosen I T, Fox E J, Jiang L, Lyu B, Li H, Watanabe K, Taniguchi T, Jung J, Shi Z, Goldhaber-Gordon D, Zhang Y and Wang F 2019 Nature 572 215–19
  * [20] Aoki M and Amawashi H 2007 Solid State Communications 142 123–7
  * [21] Nery J P, Calandra M and Mauri F 2020 Nano Letters 20 5017–23 pMID: 32525317
  * [22] Kittel C, McEuen P and McEuen P 1996 Introduction to solid state physics vol 8 (Wiley New York)
  * [23] Boehm H and Coughlin R 1964 Carbon 2 1–6
  * [24] Menezes M G, Capaz R B and Louie S G 2014 Physical Review B 89(3) 035431–5
  * [25] Lui C H, Li Z, Mak K F, Cappelluti E and Heinz T F 2011 Nature Physics 7 944–7
  * [26] Zhang F, Jung J, Fiete G A, Niu Q and MacDonald A H 2011 Physical Review Letters 106 1568014–1
  * [27] Kopnin N, Ijäs M, Harju A and Heikkilä T 2013 Physical Review B 87 140503–6
  * [28] Slizovskiy S, McCann E, Koshino M and Fal’ko V I 2019 Communications Physics 2 1–10
  * [29] Kim K S, Zhao Y, Jang H, Lee S Y, Kim J M, Kim K S, Ahn J H, Kim P, Choi J Y and Hong B H 2009 Nature 457 706–10
  * [30] Wang L, Zihlmann S, Baumgartner A, Overbeck J, Watanabe K, Taniguchi T, Makk P and Schönenberger C 2019 Nano Letters 19 4097–102
  * [31] Vincent T, Panchal V, Booth T, Power S R, Jauho A P, Antonov V and Kazakova O 2018 2D Materials 6 015022–21
  * [32] Bouhafs C, Pezzini S, Geisenhof F R, Mishra N, Mišeikis V, Niu Y, Struzzi C, Weitz R T, Zakharov A A, Forti S et al. 2021 Carbon 177 282–90
  * [33] Sanctis A D, Mehew J D, Alkhalifa S, Withers F, Craciun M F and Russo S 2018 Nano Letters 18 7919–26
  * [34] Nery J P, Calandra M and Mauri F 2021 2D Materials 8 035006–15
  * [35] Purdie D G, Pugno N M, Taniguchi T, Watanabe K, Ferrari A C and Lombardo A 2018 Nature Communications 9 5387–98
  * [36] Jiang L, Wang S, Shi Z, Jin C, Utama M I B, Zhao S, Shen Y R, Gao H J, Zhang G and Wang F 2018 Nature Nanotechnology 13 204–208
  * [37] Zhao H, Lin Y C, Yeh C H, Tian H, Chen Y C, Xie D, Yang Y, Suenaga K, Ren T L and Chiu P W 2014 ACS Nano 8 10766–73
  * [38] Huang M, Bakharev P V, Wang Z J, Biswal M, Yang Z, Jin S, Wang B, Park H J, Li Y, Qu D et al. 2020 Nature Nanotechnology 15 289–95
  * [39] Kresse G and Hafner J 1993 Physical Review B 47(1) 558–61
  * [40] Kresse G and Furthmüller J 1996 Computational Materials Science 6 15 – 50
  * [41] Kresse G and Furthmüller J 1996 Physical Review B 54(16) 11169–86
  * [42] Klimeš J, Bowler D R and Michaelides A 2009 Journal of Physics: Condensed Matter 22 022201–5


# Ultralong 100 ns Spin Relaxation Time in Graphite at Room Temperature

B. G. Márkus  Stavropoulos Center for Complex Quantum Matter, Department of
Physics and Astronomy, University of Notre Dame, Notre Dame, Indiana 46556,
USA  Institute for Solid State Physics and Optics, Wigner Research Centre for
Physics, Budapest H-1525, Hungary  Department of Physics, Institute of
Physics, Budapest University of Technology and Economics, Műegyetem rkp. 3.,
H-1111 Budapest, Hungary  M. Gmitra  Institute of Physics, Pavol Jozef Šafárik
University in Košice, Park Angelinum 9, 040 01 Košice, Slovakia  Institute of
Experimental Physics, Slovak Academy of Sciences, Watsonova 47, 04001 Košice,
Slovakia  B. Dóra  Department of Theoretical Physics, Institute of Physics and
MTA-BME Lendület Topology and Correlation Research Group Budapest University
of Technology and Economics, Műegyetem rkp. 3., H-1111 Budapest, Hungary  G.
Csősz  Department of Physics, Institute of Physics, Budapest University of
Technology and Economics, Műegyetem rkp. 3., H-1111 Budapest, Hungary  T.
Fehér  Department of Physics, Institute of Physics, Budapest University of
Technology and Economics, Műegyetem rkp. 3., H-1111 Budapest, Hungary  P.
Szirmai  Laboratory of Physics of Complex Matter, École Polytechnique Fédérale
de Lausanne, Lausanne CH-1015, Switzerland  B. Náfrádi  Laboratory of Physics
of Complex Matter, École Polytechnique Fédérale de Lausanne, Lausanne CH-1015,
Switzerland  V. Zólyomi  STFC Hartree Centre, Daresbury Laboratory, Daresbury,
Warrington WA4 4AD, United Kingdom  L. Forró  Stavropoulos Center for Complex
Quantum Matter, Department of Physics and Astronomy, University of Notre Dame,
Notre Dame, Indiana 46556, USA  Laboratory of Physics of Complex Matter, École
Polytechnique Fédérale de Lausanne, Lausanne CH-1015, Switzerland  J. Fabian
<EMAIL_ADDRESS>Department of Physics, University of Regensburg, 93040
Regensburg, Germany  F. Simon<EMAIL_ADDRESS>Department of Physics,
Institute of Physics and ELKH-BME Condensed Matter Research Group Budapest
University of Technology and Economics, Műegyetem rkp. 3., H-1111 Budapest,
Hungary  Institute for Solid State Physics and Optics, Wigner Research Centre
for Physics, Budapest H-1525, Hungary

###### Abstract

Graphite has been intensively studied, yet its electron spins dynamics remains
an unresolved problem even 70 years after the first experiments. The central
quantities, the longitudinal ($T_{1}$) and transverse ($T_{2}$) relaxation
times were postulated to be equal, mirroring standard metals, but $T_{1}$ has
never been measured for graphite. Here, based on a detailed band structure
calculation including spin-orbit coupling, we predict an unexpected behavior
of the relaxation times. We find, based on saturation ESR measurements, that
$T_{1}$ is markedly different from $T_{2}$. Spins injected with perpendicular
polarization with respect to the graphene plane have an extraordinarily long
lifetime of $100$ ns at room temperature. This is ten times more than in the
best graphene samples. The spin diffusion length across graphite planes is
thus expected to be ultralong, on the scale of $\sim 70~{}\mu$m, suggesting
that thin films of graphite — or multilayer AB graphene stacks — can be
excellent platforms for spintronics applications compatible with 2D van der
Waals technologies. Finally, we provide a qualitative account of the observed
spin relaxation based on the anisotropic spin admixture of the Bloch states in
graphite obtained from density functional theory calculations.

Spintronic devices require materials with a suitably long spin-relaxation
time, $\tau_{\text{s}}$. Carbon nanomaterials, such as graphite intercalated
compounds [1], graphene [2], fullerenes [3], and carbon nanotubes [4], have
been considered [5, 6, 7] for spintronics [8, 9, 10], as small spin-orbit
coupling (SOC) systems with low concentration of magnetic 13C nuclei which
contribute to a long $\tau_{\text{s}}$. However, experimental data and the
theory of spin-relaxation in carbon based materials face critical open
questions. Chiefly, the absolute value of $\tau_{\text{s}}$ in graphene is
debated with values ranging from $100$ ps to $12$ ns [11, 12, 13, 14, 15, 16],
and theoretical investigations suggest an extrinsic origin of the measured
short $\tau_{\text{s}}$ values [17].

Contemporary studies, in order to introduce functionality into spintronic
devices[18, 19, 20, 21, 22, 23, 24, 25, 26], focus on tailoring the SOC in
two-dimensional heterostructures with the help of proximity effect [27, 28,
29, 30]. Theory predicted a giant spin-relaxation anisotropy in graphene when
in contact with a large-SOC material [31] that was subsequently observed in
mono- and bilayer graphene [32, 33, 34, 35, 36, 37]. This is in contrast with
graphene on a SOC-free substrate having a nearly isotropic spin-relaxation
[38, 39, 40, 41]. It would be even better to have materials with an intrinsic
spin-relaxation time anisotropy, which would enable efficient control over the
spin transport and thereby boost the development of spintronic devices.

A remarkably simple example of an anisotropic carbon-based material is
graphite, which, while being one of the most extensively studied crystalline
materials, still holds several puzzles. Specifically, the spin-relaxation, its
anisotropy, and the $g$-factor are not yet understood in graphite, and this
represents a 70-year-old challenge. As early as 1953, the first spin
spectroscopic study of graphite [42, 43] used conduction electron spin
resonance (CESR). The CESR linewidth, $\Delta B$, yields directly the spin-
decoherence time: $T_{2}=\left(\gamma\Delta B\right)^{-1}$, where
$\gamma/2\pi\approx 28\,\text{GHz/T}$ is the electron gyromagnetic ratio
(which is related to the $g$-factor as $|\gamma|=g\mu_{\text{B}}/\hbar$).
Magnetic resonance is characterized by two distinct relaxation times, $T_{1}$
and $T_{2}$, which denote the relaxation of the components parallel and
perpendicular to the external magnetic field, respectively [44]. In zero
magnetic field, $T_{1}=T_{2}=\tau_{\text{s}}$ holds and the latter parameter
is measured in spin-injected transport studies.

The ESR linewidth and $g$-factor have a peculiar anisotropy in graphite.
$\Delta B$ is about twice as large for a magnetic field perpendicular to the
graphene layers (denoted as $\Delta B_{\perp}$) as compared to when the
magnetic field is in the $(a,b)$ plane (denoted as $\Delta B_{\parallel}$).
The respective $g$-factors are also strongly anisotropic: $g_{\perp}$ is
considerably shifted with respect to the free-electron value of
$g_{0}=2.0023$, and is temperature-dependent, whereas $g_{\parallel}$ is
barely shifted and is temperature independent [45, 46].

Although several explanations have been proposed [47, 45, 46, 48], no
consistent picture has emerged yet for these anomalous findings in graphite.
Nevertheless, understanding the spin-relaxation mechanism would be important
for the advancement of spin-relaxation theory in general, but especially for
spintronics applications of mono- or few-layer graphene.

We unravel the anomalous spin relaxation in graphite by studying the details
of the spin-orbit coupling, its dependence over the Fermi surface, and its
anisotropy. We find that the SOC is strongly anisotropic in graphite due to
the symmetry [47]: the (pseudospin) spin-orbit field is oriented along $z$.
This results in a significant anisotropy of the electron spin dynamics, spin-
relaxation time, and a giant $g$-factor anisotropy. We demonstrate that the
ESR data is compatible with the SOC anisotropy scenario. Remarkably, the
theory predicts an ultralong spin-relaxation time for spins that are polarized
perpendicular to the graphite $(a,b)$ plane. Indeed, saturated ESR experiments
reveal a spin-relaxation time longer than $100$ ns _at room temperature_ for
this geometry. These findings qualify graphite thin film as a strong candidate
for spintronics technology.

##  I Theoretical predictions

We performed first-principles calculations of the electronic structure of
graphite in the presence of spin-orbit coupling. Figure 1. shows the band
dispersions near the K point of the Brillouin zone (BZ). We limit our
qualitative considerations to the shown energy scale, which corresponds to the
quasiparticle energy smearing up to room temperature. The full calculated band
structure along high symmetry lines is presented in the Supporting Information
(SI).

Figure 1: Top: Calculated low-energy electronic band structure of graphite at
$K$ along the lines towards $\Gamma$ and M points. Note the band degeneracy
near the K point. Middle and bottom: Spin-mixing parameters $b_{k}^{2}$ for
the spin quantization (magnetic field $\bf B$) along $c$ and in the plane
$(a,b)$, respectively. Also indicated are the spin relaxation anisotropies for
the momenta at the vertical dashed lines, implied by the calculated
$b_{k}^{2}$ there.

In Elliott-Yafet’s theory of spin relaxation[49, 50], which is dominant in
centrosymmetric metals such as graphite, spin-relaxation is due to the mixing
of the otherwise pure spin up/down states [51, 9, 52, 53, 54, 55]. For a
selected spin quantization axis, the two degenerate Bloch states can be
characterized as up, $\ket{\uparrow}$, or down, $\ket{\downarrow}$, in the
absence of SOC. Due to spin-orbit coupling, the spin up/down states have a
(typically small) admixture of the Pauli spin down/up spinors:

$\displaystyle\ket{\widetilde{\uparrow}}_{k}=$
$\displaystyle\left[a_{k}\ket{\uparrow}+b_{k}\ket{\downarrow}\right]\mathrm{e}^{\mathrm{i}kr},$
(1) $\displaystyle\ket{\widetilde{\downarrow}}_{k}=$
$\displaystyle\left[a_{-k}^{\ast}\ket{\downarrow}+b_{-k}^{\ast}\ket{\uparrow}\right]\mathrm{e}^{\mathrm{i}kr}.$
(2)

Here $\ket{\widetilde{\uparrow}}$, and $\ket{\widetilde{\downarrow}}$ are the
Bloch states in the presence of the SOC. The spin-flip probability is
proportional to $b^{2}_{k}$.

Thanks to our high precision calculation of the dispersion relation, we could
determine $b_{k}^{2}$ near the Fermi level in Fig. 1. For the electron spin
quantized along the $c$-axis, the spin admixture probability is rather small,
less than $10^{-6}$, without a significant momentum dependence. However,
$b_{k}^{2}$ exhibits peaks at the band crossings (in fact, those are
anticrossings as SOC opens a gap of $24~{}\mu$eV such as in graphene [56])
when the spin quantization axis is in the plane. Then the spin admixture
probabilities are orders of magnitude higher than elsewhere in the Brillouin
zone. In fact, we find here so-called spin hot-spots [51], similar to what
happens in monolayer graphene at the Dirac point [57]. Fig. 1. reveals that
the dominant contribution to the spin admixture comes from the momenta along
$\textbf{K}\to\mathbf{\Gamma}$.

While calculating Fermi-level averages of $b_{k}^{2}$ is beyond the scope of
the present work (mainly due to the tiny Fermi surface of graphite and the
presence of the spin hot-spots), our calculation suggests that the spin-
relaxation rate in an in-plane magnetic field is expected to be at least an
order of magnitude faster than in a magnetic field parallel to the $c$-axis.

Figure 2: Schematics of the geometry of the different spin-relaxation
contributions for $B_{\perp}$ (upper panel) and $B_{\parallel}$ (lower panel).
Note that spins precess around the magnetic field for both orientations. The
respective spin-relaxation rates are caused by the fluctuating components
which are perpendicular to the given direction. For the $B_{\parallel}$, it
consists of two contributions due to the Larmor precession of the spins.

To relay on the spin-relaxation rate to the ESR line, we evoke the theory of
Yafet[50] developed for the case of an anisotropic SOC. For graphite, the
major features could be nicely followed in Fig. 2. He argued that both the
$T_{1}$ and $T_{2}$ relaxation times are caused by fluctuating magnetic fields
$\delta\mathbf{B}=(\delta B_{a},\delta B_{b},\delta B_{c})$ due to the SOC.
Fluctuating fields along a given direction give rise to spin relaxation of a
spin component perpendicular to them. (A more recent and rigorous derivation
of the ESR relaxation times is found in Ref. 58). Following Yafet[50], we
introduce $\delta B_{\perp}^{2}$ and $\delta B_{\parallel}^{2}$ for the
squared magnitude of the fluctuating fields along the crystalline $c$ axis,
and in the $(a,b)$ plane, respectively. The results for the corresponding
spin-relaxation rates are:

$\displaystyle\left(T_{1}\right)^{-1}_{\perp}$ $\displaystyle\propto 2\delta
B_{\parallel}^{2}$ (3) $\displaystyle\left(T_{2}\right)^{-1}_{\perp}$
$\displaystyle\propto\delta B_{\parallel}^{2}+\delta B_{\perp}^{2}$ (4)
$\displaystyle\left(T_{1}\right)^{-1}_{\parallel}$ $\displaystyle\propto\delta
B_{\parallel}^{2}+\delta B_{\perp}^{2}$ (5)
$\displaystyle\left(T_{2}\right)^{-1}_{\parallel}$
$\displaystyle\propto\frac{3}{2}\delta B_{\parallel}^{2}+\frac{1}{2}\delta
B_{\perp}^{2}.$ (6)

This is illustrated in Fig. 2. E.g., for $B_{\perp}$, spins precess around the
$c$ axis and $T_{1}$ is caused entirely by $\delta B_{\parallel}$ (Eq. (3)),
however $T_{2}$ is caused by both $\delta B_{\parallel}$ and $\delta
B_{\perp}$ fluctuating fields as these SOC fields are perpendicular to spin
component which is in the plane (Eq. (4)). For the $B_{\parallel}$
orientation, the spins precess in the $a-c$ plane thus $T_{1}$ is caused by
both the $\delta B_{\parallel}$ and $\delta B_{\perp}$ which also yields
$T_{2,\perp}=T_{1,\parallel}$ (Eqs. (4) and (5)). Yafet discussed the case of
the _extreme uniaxial anisotropy_ , i.e., when $\delta B_{\parallel}^{2}=0$
while $\delta B_{\perp}^{2}$ is finite, which is in fact, due to symmetry
considerations, the theoretical prediction for graphene (Ref. 59). We suppose
that this is the case for graphite, as well, which we check experimentally. It
is interesting to notice, that somewhat counterintuitively, the extreme
anisotropy corresponds to just a factor of $2$ anisotropy in the ESR linewidth
as $\left(T_{2}\right)^{-1}_{\perp}/\left(T_{2}\right)^{-1}_{\parallel}=\Delta
B_{\perp}/\Delta B_{\parallel}=2$ (see Eqs. (4) and (6)).

##  II Results and Discussion

###  II.1 Anisotropy of the ESR linewidth and $g$-factor

Figure 3: Temperature-dependent ESR data at $9.4$ GHz ($\sim 0.3$ T) in
graphite for both crystallographic orientations. Upper panel: the ESR
linewidth with different symbols for the as-obtained and annealed samples. For
$B_{\perp}$, the data is shown only above $30$ K. Lower panel: linewidth
anisotropy factor, i.e. the ratio of the linewidth for $B_{\perp}$ and
$B_{\parallel}$ configurations. The dashed line is the constant $2$.

Fig. 3 shows the temperature dependence of the ESR linewidth for the two major
magnetic field orientations in highly-oriented pyrolytic graphite (HOPG)
samples. Comparing the data with previous measurements[43, 60, 61, 62, 46], we
find good agreement for the $50-300$ K temperature range. Due to the finite
penetration depth of microwaves, the ESR lineshape is asymmetric. This is
taken into account by fitting the ESR spectra to the so-called Dysonian lines
[63] that are well approximated by a mixture of absorption and dispersion
components of a Lorentzian curve [64]. Additional data including typical ESR
spectra, the $g$-factor, the temperature-dependent ESR intensity, angular
dependence of the ESR spectra, and linewidth are provided in the SI.

Above $50$ K, the linewidth increases with decreasing temperature for both
orientations. Linewidth data for $B_{\perp}$ are not shown below $30$ K where
the significant line broadening and the relatively weak signal intensity
prevent a reliable analysis. Annealing at $300~{}^{\circ}\text{C}$ in a
dynamic vacuum affects the ESR linewidth only below $50$ K. This is plausible
since due to its semi-metallic character and the small Fermi surface even a
few ppm of paramagnetic impurities can dominate the signal at low temperatures
due to their Curie spin susceptibility[65]. Higher annealing temperatures (up
to $600~{}^{\circ}\text{C}$) did not influence further the linewidth.

From a calibrated measurement of the spin-susceptibility, we obtained that the
Curie contribution to the ESR signal is due to 5(1) ppm of $S=1/2$ spins. We
show in Fig. S2 that the susceptibility contribution from the localized and
delocalized spins is equal at about 15 K. Given that any influence of the
localized species on the $g$-factor and the linewidth is weighted with the
spin-susceptibility [65], we conclude that above $50-100$ K (especially at the
technologically important room temperature region), the observed ESR linewidth
is the intrinsic property of graphite. Although the present work focuses on
the anisotropy of the ESR linewidth, we also mention that the conventional
Elliott-Yafet’s theory predicts an opposite temperature dependence, which may
be related to the near band degeneracy around the graphite Fermi surface. Its
full description requires additional work though.

The temperature-dependent $g$-factor data (see Fig. S5 of SI) confirms the
earlier experimental data[60, 61, 62, 46] and the theoretical prediction. The
measured factor of $2$ anisotropy in the ESR linewidth of graphite (see the
bottom image of Fig. 3) matches _exactly_ the theoretical prediction of a
vanishing SOC for spins polarized in the $(a,b)$ plane and a finite SOC for
spins polarized along the $c$-axis. Strictly speaking, the linewidth
anisotropy ratio is $2$ only between $100$ and $300$ K, while it deviates
somewhat upwards below $100$ K and downwards above $300$ K, which calls for
more detailed modeling.

###  II.2 The ultralong spin-relaxation time, $T_{1}$

An important consequence of the absence of spin-orbit fields in the graphene
planes is the predicted ultralong spin-lattice relaxation time,
$T_{1,{\perp}}$, according to Eq. (3), whereas the other three relaxation
times remain finite, below $50$ ns. To experimentally verify this prediction,
we performed saturation ESR measurements (see Fig. S7 of SI). Although the
transversal spin-relaxation time, $T_{2}$ is accessible directly from the
linewidth, $T_{1}$ does not directly affect the lineshape in conventional ESR
studies. In principle, both relaxation times can be measured with a spin-echo
technique but it is limited to relaxation times when both $T_{1}$ and $T_{2}$
are longer than $\sim 1~{}\mu$s, which is not the case herein. However,
saturation ESR experiments [44] yield relaxation times down to $\sim 10$ ns.
The method is based on monitoring the variation of the ESR signal intensity
and the linewidth as a function of the irradiating microwave power. A
characteristic drop in the signal intensity, accompanied by a line broadening
allows determining $T_{1}$. Further technical details are given in the SI
together with the raw saturation ESR data.

Figure 4: The spin-lattice relaxation time, $T_{1}$, for the two orientations
of the magnetic field (red triangles: $B_{\parallel}$, blue dots: $B_{\perp}$)
in graphite as a function of temperature. Note the about $10$ times longer
$T_{1}$ when $B_{\perp}$ as compared to $B_{\parallel}$. The inset depicts the
expected experimental situation in a spin transport experiment. Dashed line
denotes $T_{2}$ when $B_{\perp}$, which also confirms the case of extreme
anisotropy as it aligns well with the $T_{1}$ values for $B_{\parallel}$, as
predicted by Eqs. (4) and (5).

The $T_{1}$ results for both orientations of the magnetic field are shown in
Fig. 4. The data is given only above $300$ K for $B_{\perp}$ and above $370$ K
for $B_{\parallel}$ as below these temperatures the shortening of the
respective $T_{2}$ prevents the observation of the saturation effect.
Nevertheless, having in mind the spintronics application, this is the relevant
temperature range. Red triangles and the dashed line in Fig. 4 confirm that
the $T_{1,\parallel}~{}=~{}T_{2,\perp}$ prediction from Eqs. (4) and (5) is
indeed satisfied, which provides further proof that the extreme SOC anisotropy
is realized in graphite.

The most striking experimental observation is the presence of $T_{1,\perp}$
relaxation times beyond $100$ ns and the approximate factor $10$ anisotropy of
$T_{1}$. Although it is observed in an ESR experiment, i.e., in the presence
of an external magnetic field, this result can be readily extended to the case
of zero magnetic fields: it implies that $\tau_{\text{s},\perp}$ in graphite
would be as long as $100$ ns when electrons are injected with a spin
perpendicular to the graphene planes and it would be about a factor $10$ times
shorter when the spins lie in the graphene planes. This is also depicted in
Fig. 4. The giant anisotropy may find a number of applications, e.g., it could
be exploited to control the spin-relaxation times in spintronic devices.

It is known [1] that electrons travel mainly in the graphene planes in
graphite and interlayer transport is diffusion limited and it leads to a
conduction anisotropy $\sigma_{\parallel}/\sigma_{\perp}$ beyond
$10^{3}-10^{4}$. This allows to approximate the spin diffusion in graphite
while considering that electrons reside on a given graphene layer. The large
value of $\tau_{\text{s},\perp}$ leads to an ultra-long spin-diffusion length
in graphite:
$\delta_{\text{s}}=\frac{v_{\text{F}}}{\sqrt{2}}\sqrt{\tau\tau_{\text{s}}}$
(the factor $2$ appears from the two-dimensional diffusion equation). It gives
$\delta_{\text{s}}\approx 70~{}\mu\text{m}$ with typical values of
$v_{\text{F}}=10^{6}$ m$/$s and $\tau=10^{-13}$ s. This is already a
macroscopic length scale, which could bring spintronic devices closer to
reality.

##  III Summary

In conclusion, we unraveled the anomalous dynamics of itinerant electron spins
in graphite. We first studied the details of the band-structure-related spin-
orbit coupling which predicted a hitherto hidden, symmetry-related extreme
anisotropy of the spin-orbit coupling. This anisotropy in fact explains the
known anomalous anisotropic properties of graphite: a strongly anisotropic
$g$-factor and the $1:2$ anisotropy of the ESR linewidth. We recognized that
the latter should in principle be accompanied by a giant anisotropy of the
spin-lattice relaxation time. This prediction was examined with saturation ESR
studies and we observe ultralong ($T_{1}$ in excess of $100$ ns) spin-
relaxation times for spins aligned perpendicular to the graphene planes. When
extrapolated to the zero-field limit, this predicts similarly long spin-
relaxation times in spin-transport studies and a macroscopic spin-diffusion
length.

## Methods

#### Experiment

We studied high-quality HOPG (highly-oriented pyrolitic graphite) from
Structure Probe Inc. (SPI Grade I) using electron spin resonance (ESR). The
HOPG had a mosaicity of $0.4^{\circ}\pm 0.1^{\circ}$. We studied X-band
($0.33$ T, $9.4$ GHz) in a commercial ESR spectrometer (Bruker Elexsys E500)
in the $4-673$ K temperature range. The spectrometer is equipped with a
goniometer which allows reliable and reproducible sample rotations. We used
HOPG disk samples of $3$ mm diameter and $70~{}\mu\text{m}$ thickness.

The samples were sealed in quartz ampules under $20$ mbar He for the ESR
measurement. Annealing of some of the samples was performed in a high dynamic
vacuum in a furnace up to $300-600~{}^{\circ}$C. The goals in this study are
the accurate measurement of the ESR linewidth, and temperature-dependent ESR
intensity (penetration effects for a bulky sample prevent the determination of
the absolute spin-susceptibility). We also monitored the $g$-factor as a
function of temperature to enable comparison with existing literature data. We
employed a low magnetic field modulation to avoid signal distortion, we also
used a low ($150~{}\mu\text{W}$) microwave power to avoid saturation of the
normal data and we monitored the change in the cavity resonance frequency. We
quote the widths of derivative Lorentzian curves (or half width-half maximum
data, HWHM) fitted to the experimental data which is related to the peak-to-
peak linewidth as: $\Delta B_{\text{HWHM}}=\Delta B_{\text{pp}}/\sqrt{3}$.
Saturated ESR experiments were performed in a commercial microwave cavity
(Bruker ER 4122 SHQ, Super High $Q$ Resonator) with an unloaded $Q_{0}=7500$
and with the sample $Q_{\text{L}}=5500$. This cavity produces AC magnetic
fields of $B_{1}=0.2~{}\text{mT}\sqrt{\frac{pQ_{\text{L}}}{Q_{0}}}$, where $p$
is the microwave power in Watts and we used microwave powers up to $0.2$ W.

#### Theory

The electronic band structure of graphite in the presence of spin-orbit
coupling is determined by the full potential linearized augmented plane waves
(LAPW) method, based on density functional theory implemented in Wien2k [66].
For exchange-correlation effects, the generalized gradient approximation was
utilized [67]. In our three-dimensional calculation the graphene sheets of
lattice constant $a=1.42\sqrt{3}$ Å are separated by the distance of $c=3.35$
Å. Integration in the reciprocal space was performed by the modified Blöchl
tetrahedron scheme, taking the mesh of $33\times 33$ $k$-points in the
irreducible Brillouin zone (BZ) wedge. As the plane-wave cut-off, we took
$9.87$ Å-1. The $1$s core states were treated fully relativistically by
solving the Dirac equation, while spin-orbit coupling for the valence
electrons was treated within the muffin-tin radius of $1.34$ a.u. by the
second variational method [68].

The algorithm and the code used to calculate the $b_{k}^{2}$ values were
proven to be correct for various other materials previously [69]. These
include WS2[70], WSe2, MoSe2[71], various heterostructures [72, 73], etc. We
are confident that it provides adequate results for the case of graphite as
well. However, averaging the spin admixture, $b^{2}$, over the tiny Fermi-
surface pockets is, at the moment, not feasible. The main reason is that the
admixture varies strongly (over several orders of magnitude) close to the band
degeneracies (spin hot spots), requiring very dense sampling of the Fermi
surface. At the moment, this is beyond the computational capabilities of our
available computing infrastructure.

## Data availability

The data needed to evaluate and reproduce the conclusions are present in the
paper and the Supporting Information (Online Content). Additional data related
to this paper are available from the corresponding author upon request.

## Online Content

Online content includes supporting information.

## References

  * [1] Dresselhaus, M. S. & Dresselhaus, G.  Intercalation compounds of graphite.  _Adv. Phys._ 51, 1–186 (2002).
  * [2] Novoselov, K. S. _et al._ Electric field effect in atomically thin carbon films.  _Science_ 306, 666–669 (2004).
  * [3] Kroto, H. W., Heath, J. R., O’Brien, S. C., Curl, R. F. & Smalley, R. E.  C60: Buckminsterfullerene.  _Nature_ 318, 162–163 (1985).
  * [4] Iijima, S.  Helical microtubules of graphitic carbon.  _Nature_ 354, 56–58 (1991).
  * [5] Han, R. K., W. Kawakami, Gmitra, M. & Fabian, J.  Graphene spintronics.  _Nature Nanotechnology_ 9, 794–807 (2014).
  * [6] Roche, S. _et al._ Graphene spintronics: the European Flagship perspective.  _2D Materials_ 2, 030202 (2015).
  * [7] Avsar, A. _et al._ Colloquium: Spintronics in graphene and other two-dimensional materials.  _Rev. Mod. Phys._ 92, 021003 (2020).
  * [8] Wolf, S. A. _et al._ Spintronics: A Spin-Based Electronics Vision for the Future.  _Science_ 294, 1488–1495 (2001).
  * [9] Žutić, I., Fabian, J. & Sarma, S. D.  Spintronics: Fundamentals and applications.  _Rev. Mod. Phys._ 76, 323–410 (2004).
  * [10] Wu, M. W., Jiang, J. H. & Weng, M. Q.  Spin dynamics in semiconductors.  _Phys. Rep._ 493, 61–236 (2010).
  * [11] Tombros, N., Józsa, C., Popinciuc, M., Jonkman, H. T. & van Wees, B. J.  Electronic spin transport and spin precession in single graphene layers at room temperature.  _Nature_ 448, 571–574 (2007).
  * [12] Han, W. & Kawakami, R. K.  Spin Relaxation in Single-Layer and Bilayer Graphene.  _Phys. Rev. Lett._ 107, 047207 (2011).
  * [13] Yang, T.-Y. _et al._ Observation of Long Spin-Relaxation Times in Bilayer Graphene at Room Temperature.  _Phys. Rev. Lett._ 107, 047206 (2011).
  * [14] Roche, S. & Valenzuela, S. O.  Graphene spintronics: puzzling controversies and challenges for spin manipulation.  _Journal of Physics D: Applied Physics_ 47, 094011 (2014).
  * [15] Venkata Kamalakar, M., Groenveld, C., A., D. & Dash, S. P.  Long distance spin communication in chemical vapour deposited graphene.  _Nature Communications_ 6, 6766 (2015).
  * [16] Drögeler, M. _et al._ Spin Lifetimes Exceeding 12 ns in Graphene Nonlocal Spin Valve Devices.  _Nano Letters_ 16, 3533–3539 (2016).
  * [17] Kochan, D., Gmitra, M. & Fabian, J.  Spin Relaxation Mechanism in Graphene: Resonant Scattering by Magnetic Impurities.  _Phys. Rev. Lett._ 112, 116602 (2014).
  * [18] Avsar, A. _et al._ Spin–orbit proximity effect in graphene.  _Nature Communications_ 5, 4875 (2014).
  * [19] Wang, Z. _et al._ Strong interface-induced spin–orbit interaction in graphene on WS2.  _Nature Communications_ 6, 8339 (2015).
  * [20] Wang, Z. _et al._ Origin and Magnitude of ’Designer’ Spin-Orbit Interaction in Graphene on Semiconducting Transition Metal Dichalcogenides.  _Phys. Rev. X_ 6, 041020 (2016).
  * [21] Yang, B. _et al._ Tunable spin–orbit coupling and symmetry-protected edge states in graphene/WS2.  _2D Materials_ 3, 031012 (2016).
  * [22] Yang, B. _et al._ Strong electron-hole symmetric Rashba spin-orbit coupling in graphene/monolayer transition metal dichalcogenide heterostructures.  _Phys. Rev. B_ 96, 041409 (2017).
  * [23] Torres, W. S. _et al._ Spin precession and spin Hall effect in monolayer graphene/Pt nanostructures.  _2D Materials_ 4, 041008 (2017).
  * [24] Offidani, M., Milletarì, M., Raimondi, R. & Ferreira, A.  Optimal Charge-to-Spin Conversion in Graphene on Transition-Metal Dichalcogenides.  _Phys. Rev. Lett._ 119, 196801 (2017).
  * [25] Dankert, A. & Dash, S. P.  Electrical gate control of spin current in van der Waals heterostructures at room temperature.  _Nature Communications_ 8, 16093 (2017).
  * [26] Benítez, L. A. _et al._ Tunable room-temperature spin galvanic and spin Hall effects in van der Waals heterostructures.  _Nat. Mat._ 19, 170–175 (2020).
  * [27] Gmitra, M. & Fabian, J.  Graphene on transition-metal dichalcogenides: A platform for proximity spin-orbit physics and optospintronics.  _Phys. Rev. B_ 92, 155403 (2015).
  * [28] Gmitra, M. & Fabian, J.  Proximity Effects in Bilayer Graphene on Monolayer ${\mathrm{WSe}}_{2}$: Field-Effect Spin Valley Locking, Spin-Orbit Valve, and Spin Transistor.  _Phys. Rev. Lett._ 119, 146401 (2017).
  * [29] Ẑutić, I., A., M.-A., B., S., Dery, H. & Belashchenko, K.  Proximitized materials”.  _Materials Today_ 22, 85 – 107 (2019).
  * [30] Högl, P. _et al._ Quantum Anomalous Hall Effects in Graphene from Proximity-Induced Uniform and Staggered Spin-Orbit and Exchange Coupling.  _Phys. Rev. Lett._ 124, 136403 (2020).
  * [31] Cummings, A. W., Garcia, J. H., Fabian, J. & Roche, S.  Giant Spin Lifetime Anisotropy in Graphene Induced by Proximity Effects.  _Phys. Rev. Lett._ 119, 206601 (2017).
  * [32] Benítez, L. A. _et al._ Strongly anisotropic spin relaxation in graphene–transition metal dichalcogenide heterostructures at room temperature.  _Nat. Phys._ 14, 303–308 (2018).
  * [33] Leutenantsmeyer, J. C., Ingla-Aynés, J., Fabian, J. & van Wees, B. J.  Observation of Spin-Valley-Coupling-Induced Large Spin-Lifetime Anisotropy in Bilayer Graphene.  _Phys. Rev. Lett._ 121, 127702 (2018).
  * [34] Zihlmann, S. _et al._ Large spin relaxation anisotropy and valley-Zeeman spin-orbit coupling in ${\mathrm{WSe}}_{2}$/graphene/$h$-BN heterostructures.  _Phys. Rev. B_ 97, 075434 (2018).
  * [35] Xu, J., Zhu, T., Luo, Y. K., Lu, Y.-M. & Kawakami, R. K.  Strong and Tunable Spin-Lifetime Anisotropy in Dual-Gated Bilayer Graphene.  _Phys. Rev. Lett._ 121, 127703 (2018).
  * [36] Wakamura, T. _et al._ Strong Anisotropic Spin-Orbit Interaction Induced in Graphene by Monolayer ${\mathrm{WS}}_{2}$.  _Phys. Rev. Lett._ 120, 106802 (2018).
  * [37] Omar, S., Madhushankar, B. N. & van Wees, B. J.  Large spin-relaxation anisotropy in bilayer-graphene/${\mathrm{WS}}_{2}$ heterostructures.  _Phys. Rev. B_ 100, 155415 (2019).
  * [38] Tombros, N. _et al._ Anisotropic Spin Relaxation in Graphene.  _Phys. Rev. Lett._ 101, 046601 (2008).
  * [39] Raes, B. _et al._ Determination of the spin-lifetime anisotropy in graphene using oblique spin precession.  _Nature Communications_ 7, 11444 (2016).
  * [40] Zhu, T. & Kawakami, R. K.  Modeling the oblique spin precession in lateral spin valves for accurate determination of the spin lifetime anisotropy: Effect of finite contact resistance and channel length.  _Phys. Rev. B_ 97, 144413 (2018).
  * [41] Ringer, S. _et al._ Measuring anisotropic spin relaxation in graphene.  _Phys. Rev. B_ 97, 205439 (2018).
  * [42] Castle, J. G.  Paramagnetic Resonance Absorption in Graphite.  _Phys. Rev._ 92, 1063–1063 (1953).
  * [43] Hennig, G. R., Smaller, B. & Yasaitis, E. L.  Paramagnetic Resonance Absorption in Graphite.  _Phys. Rev._ 95, 1088–1089 (1954).
  * [44] Slichter, C. P.  _Principles of Magnetic Resonance_ (Spinger-Verlag, New York, 1989), 3rd ed. 1996 edn.
  * [45] Sercheli, M. S., Kopelevich, Y., Ricardo da Silva, R., Torres, J. H. S. & Rettori, C.  Evidence for internal field in graphite: a conduction electron spin-resonance study.  _Solid State Communications_ 121, 579–583 (2002).
  * [46] Huber, D. L., Urbano, R. R., Sercheli, M. S. & Rettori, C.  Fluctuating field model for conduction electron spin resonance in graphite.  _Phys. Rev. B_ 70, 125417 (2004).
  * [47] McClure, J. W. & Yafet, Y.  Theory of the $g$-factor of the current carriers in graphite single crystals.  _Proc. 5th Conference on Carbon, Pergamon Press/The Macmillan Company, New York_ 22–28 (1962).
  * [48] Restrepo, O. D. & Windl, W.  Full First-Principles Theory of Spin Relaxation in Group-IV Materials.  _Phys. Rev. Lett._ 109, 166604 (2012).
  * [49] Elliott, R. J.  Theory of the Effect of Spin-Orbit Coupling on Magnetic Resonance in Some Semiconductors.  _Phys. Rev._ 96, 266–279 (1954).
  * [50] Yafet, Y.  $g$-Factors and Spin-Lattice Relaxation of Conduction Electrons.  _Solid State Physics_ 14, 1–98 (1963).
  * [51] Fabian, J. & Das Sarma, S.  Spin Relaxation of Conduction Electrons in Polyvalent Metals: Theory and a Realistic Calculation.  _Phys. Rev. Lett._ 81, 5624–5627 (1998).
  * [52] Kurpas, M., Gmitra, M. & Fabian, J.  Spin-orbit coupling and spin relaxation in phosphorene: Intrinsic versus extrinsic effects.  _Phys. Rev. B_ 94, 155423 (2016).
  * [53] Zimmermann, B. _et al._ Anisotropy of Spin Relaxation in Metals.  _Phys. Rev. Lett._ 109, 236603 (2012).
  * [54] Long, N. H. _et al._ Spin-flip hot spots in ultrathin films of monovalent metals: Enhancement and anisotropy of the Elliott-Yafet parameter.  _Phys. Rev. B_ 88, 144408 (2013).
  * [55] Zimmermann, B. _et al._ Fermi surfaces, spin-mixing parameter, and colossal anisotropy of spin relaxation in transition metals from ab initio theory.  _Phys. Rev. B_ 93, 144403 (2016).
  * [56] Gmitra, M., Konschuh, S., Ertler, C., Ambrosch-Draxl, C. & Fabian, J.  Band-structure topologies of graphene: Spin-orbit coupling effects from first principles.  _Phys. Rev. B_ 80, 235431 (2009).
  * [57] Kurpas, M., Faria Junior, P. E., Gmitra, M. & Fabian, J.  Spin-orbit coupling in elemental two-dimensional materials.  _Phys. Rev. B_ 100, 125422 (2019).
  * [58] Fabian, J., Matos-Abiaguea, A., Ertlera, C., Stano, P. & Zutic, I.  Semiconductor Spintronics.  _Acta Physica Slovaca_ 57, 565–907 (2007).
  * [59] Kane, C. L. & Mele, E. J.  Quantum Spin Hall Effect in Graphene.  _Phys. Rev. Lett._ 95, 226801 (2005).
  * [60] Wagoner, G.  Spin Resonance of Charge Carriers in Graphite.  _Phys. Rev._ 118, 647–653 (1960).
  * [61] Singer, L. S. & Wagoner, G.  Electron Spin Resonance in Polycrystalline Graphite.  _The Journal of Chemical Physics_ 37, 1812–1817 (1962).
  * [62] Matsubara, K., Tsuzuku, T. & Sugihara, K.  Electron spin resonance in graphite.  _Phys. Rev. B_ 44, 11845–11851 (1991).
  * [63] Dyson, F. J.  Electron spin resonance absorption in metals II. Theory of electron diffusion and the skin effect.  _Phys. Rev._ 98, 349–359 (1955).
  * [64] Walmsley, L.  Translating Conduction–Electron Spin–Resonance Lines into Lorentzian Lines.  _Journal of Magnetic Resonance, Series A_ 122, 209–213 (1996).
  * [65] Barnes, S. E.  Theory of electron spin resonance of magnetic ions in metals.  _Adv. Phys._ 30, 801–938 (1981).
  * [66] Blaha, P. _et al._ WIEN2k: An APW+lo program for calculating the properties of solids.  _J. Chem. Phys._ 152, 074101 (2020).
  * [67] Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized Gradient Approximation Made Simple.  _Phys. Rev. Lett._ 77, 3865 (1996).
  * [68] Singh, D. L. & Nordstrom, L.  _Planewaves, Pseudopotentials, and LAPW method_ (Spinger-Verlag, 2006).
  * [69] Junior, P. E. F. _et al._ First-principles insights into the spin-valley physics of strained transition metal dichalcogenides monolayers.  _New J. Phys._ 24, 083004 (2022).
  * [70] Blundo, E. _et al._ Strain-Induced Exciton Hybridization in ${\mathrm{WS}}_{2}$ Monolayers Unveiled by Zeeman-Splitting Measurements.  _Phys. Rev. Lett._ 129, 067402 (2022).
  * [71] Raiber, S. _et al._ Ultrafast pseudospin quantum beats in multilayer WSe2 and MoSe2.  _Nat. Commun._ 13, 4997 (2022).
  * [72] Zollner, K. & Fabian, J.  Engineering Proximity Exchange by Twisting: Reversal of Ferromagnetic and Emergence of Antiferromagnetic Dirac Bands in $\text{Graphene/}{\mathrm{Cr}}_{2}{\mathrm{Ge}}_{2}{\mathrm{Te}}_{6}$.  _Phys. Rev. Lett._ 128, 106401 (2022).
  * [73] Zollner, K., Faria Junior, P. E. & Fabian, J.  Strong manipulation of the valley splitting upon twisting and gating in ${\mathrm{MoSe}}_{2}/{\mathrm{CrI}}_{3}$ and ${\mathrm{WSe}}_{2}/{\mathrm{CrI}}_{3}$ van der Waals heterostructures.  _Phys. Rev. B_ 107, 035112 (2023).

## Acknowledgements

Work supported by the National Research, Development and Innovation Office of
Hungary (NKFIH) Grants Nr. K137852, K142179, 2022-2.1.1-NL-2022-00004, and
2019-2.1.7-ERA-NET-2021-00028. The Swiss National Science Foundation (Grant
No. 200021 144419), DFG (German Research Foundation), SFB 1277 (Project No.
314695032), SPP 2244 (Project No. 443416183), European Union Horizon 2020
Research and Innovation Program under contract number 881603 (Graphene
Flagship), and the project FLAG ERA JTC 2021 2DSOTECH are acknowledged. M.G.
acknowledges financial support from Slovak Research and Development Agency
provided under Contract No. APVV-SK-CZ-RD-21-0114 and by the Ministry of
Education, Science, Research and Sport of the Slovak Republic provided under
Grant No. VEGA 1/0105/20 and Slovak Academy of Sciences project IMPULZ
IM-2021-42.

## Author contributions

BGM, PSz, BN, and TF performed the ESR studies under the supervision of LF.
BGM performed the high temperature saturated ESR studies. MG, VZ, and JF
performed the band structure calculations and the analysis of the SOC. BD and
FS developed the model for the temperature dependent ESR linewidth and
$g$-factor. GCs performed numerical simulations of the ESR lineshape. FS and
JF outlined the overall explanation for the spin-relaxation properties. All
authors contributed to the writing of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Supporting information is available for this paper.