# Unfaulting mechanisms of Frank loops in fluorite oxides

Miaomiao Jin  Corresponding author: M. Jin<EMAIL_ADDRESS>Department of
Nuclear Engineering, The Pennsylvania State University, 218 Hallowell
Building, University Park, 16802 PA, USA  Jilang Miao  Department of Nuclear
Engineering, The Pennsylvania State University, 218 Hallowell Building,
University Park, 16802 PA, USA  Yongfeng Zhang  Department of Engineering
Physics, University of Wisconsin-Madison, 1500 Engineering Dr, Madison, WI
53706, USA  Marat Khafizov  Department of Mechanical and Aerospace
Engineering, The Ohio State University, 201 W 19th Ave, Columbus, OH 43210,
USA  Kaustubh K. Bawane  Idaho National Laboratory, 2525 Fremont Ave, Idaho
Falls, ID 83402, USA  Boopathy Kombaiah  Idaho National Laboratory, 2525
Fremont Ave, Idaho Falls, ID 83402, USA  David H. Hurley  Idaho National
Laboratory, 2525 Fremont Ave, Idaho Falls, ID 83402, USA

###### Abstract

Unfaulting of Frank loops in irradiated fluoride oxides are of significance to
microstructural evolution. However, the mechanisms have not been directly
observed. To this end, we utilize molecular dynamics to reveal the atomistic
details related to the unfaulting process of interstitial Frank loop in ThO2,
which involve nucleation of single or multiple Shockley partial pairs at the
loop circumference. The unfaulting is achieved via a synchronous shear of the
partial pairs to remove the extrinsic stacking fault in the cation sublattice
and the intrinsic stacking fault in the anion sublattice. The strong oxygen
motion at the dislocation core may reduce the activation barriers of
dislocation nucleation and migration. These findings provide a fundamental
understanding of the transformation of faulted loops in irradiated ThO2, and
could be transferable to other fluorite systems.

Keywords

Dislocation loop, unfaulting, fluorite structure

Declaration

No conflict of interest

###### pacs:

Valid PACS appear here

Aggregation of radiation-induced defects can form dislocation loops, which are
commonly recognized under transmission electron microscopes (TEM) in
irradiated crystalline specimens chen2020effects ; kubota2005transition . In
fluorite-structured compounds, such as UO2, ThO2, CeO2, and ZrO2, the
character of dislocation loops identified in irradiated samples depends on the
irradiation conditions. In UO2 irradiated with heavy ions, only perfect loops
($\mathbf{b}=a_{0}/2\langle 110\rangle$, where $a_{0}$ is the lattice
parameter) are reported, while neutron irradiation leads to both perfect and
faulted loops ($\mathbf{b}=1/3\langle 111\rangle)$ loops. In Kr-irradiated
ThO2, both perfect and faulted loops are observed he2022dislocation . In Kr-
and Xe-irradiated CeO2, only faulted loops are confirmed
chen2013characterization . In ZrO2, faulted loops were identified under high-
voltage electron microscope baufeld1993situ . These loops are recognized as
interstitial loops, and the formation is believed to be mainly due to the
diffusion and clustering of radiation-induced cation defects
chen2013characterization ; he2022dislocation . In some other cases with
electron irradiation, faulted loops consisting of only oxygen interstitials
were also proposed in CeO2 yasunaga2008electron , however, such a high
unbalance of local charge may not be energetically stable ryazanov2002growth .

Elastic theory indicates that the perfect loop is more energetically favorable
than the faulted loop at large sizes due to the energy penalty from the
stacking fault hull2001introduction , causing faulted loops thermodynamically
unstable as size increases. As the perfect loop resultant from loop unfaulting
is mobile to interact with other defects (e.g., dislocation networks, grain
boundaries, voids, etc), the unfaulting process is of significance to
understanding the microstructure evolution. In ZrO2 at 200 oC, faulted loops
show fast growth and instantaneously transform into a perfect loop under high-
voltage electron microscope baufeld1993situ . In Kr-irradiated ThO2, the
transformation from faulted loops to perfect loops is implied in the loop
statistics: i) loop size tends to saturate and ii) the loop density tends to
decrease with increasing dose. Although various indirect results suggest the
unfaulting, to the authors’ knowledge, the unfaulting process in fluorite
oxides has not yet been directly resolved in experiments and simulations.

As the microscopic details remain largely challenging in experiments,
molecular dynamics (MD) simulations have been frequently used to examine the
loop unfaulting process in various lattice structures kadoyoshi2007molecular ;
chen2021comprehensive ; dai2022atomistic ; chen2022combined . Based on the
understanding of face-centered cubic (FCC) metals, which are pertinent to
fluorite structures, to remove the stacking fault enclosed by the loop, there
are various mechanisms involving moving Shockley partials
($\mathbf{b}=a/6<112>$) to sweep the fault chen2022combined . Removal of
intrinsic stacking fault (SF) can be realized via a single sweep of Shockley
partial, while for extrinsic SF, double sweeps of Shockley partials are
required. The Shockley partials may nucleate at the loop line or within the
loop. The resultant dislocation structure depends on the loop habit plane,
stacking fault energy, and external stress chen2022combined . Although the
understanding of FCC systems is largely established, one key issue to resolve
is whether such understanding is transferable to the fluorite oxides, where
oxygen atoms may mediate the unfaulting process.

To pursue a mechanistic understanding, we utilize large-scale MD simulations
to probe the details of loop unfaulting in ThO2, which is expected to apply to
other fluorite systems. Here we use ThO2 as a demonstration to alleviate the
concern that unpaired $f$ electrons as in UO2 and CeO2 may cause serious
inaccuracy in the interatomic potential with the loop structures. The
classical MD method as implemented in LAMMPS plimpton2007lammps  is used with
an embedded-atom method potential for ThO2 by Cooper et al. cooper2014many .
This potential takes a fixed non-formal charge model for both Th (+2.2208) and
O (-1.1104) atoms, and can well reproduce mechanical properties and defect
energies as compared to first-principles calculations and experimental
measurements. Periodic boundary conditions, and Nose-Hoover thermostat
nose1984unified ; hoover1985canonical  (if applicable) are utilized. All the
constructed systems are charge neutral to avoid numerical issues resulting
from periodic boundary conditions.

To be consistent with neutron/ion-irradiated experiments, we only focus on
loops of interstitial nature. To enable a strong driving force to initiate the
loop transformation within the MD timescale, a large faulted loop (20 nm in
diameter) is created, with a total of 974,685 atoms and a box size of around
40 nm $\times$ 40 nm $\times$ 10 nm (FIG. 1a). This loop size guarantees a
higher formation energy than that of a perfect loop, and such size is scarcely
observed in ThO2 he2022dislocation . In fact, transformation may start to
occur in nanometer-sized loops, as suggested in UO2 chartier2016early ;
le2016empirical . The loop is created by inserting three layers, i.e., O-Th-O
111 layers; such stacking has been suggested to be energetically favorable in
other fluorite oxides, e.g., UO2 le2016empirical ; aidhy2011comparison , CeO2
aidhy2011comparison ; chen2013characterization  and ZrO2 baufeld1993situ . The
stacking of layers is carefully controlled to ensure correct stacking (FIG.
1b). Although hexagonal-shaped loops are frequently observed in FCC metals,
the loops identified in fluoride oxides commonly appear in circular shape
baufeld1993situ ; he2013situ , which is hence assumed in the current study.
Atom visualization is accomplished with the OVITO package
stukowski2009visualization .

Figure 1: (a) Faulted dislocation loop configuration in ThO2. (b) Enlarged
local atom arrangement, showing the local stacking of Th and O layers. (c)
Binding energies of individual oxygen atoms in the loop (loop Th atoms are not
shown).

In the current calculations, loop stoichiometry (i.e., each inserted layer has
the exact same number of atoms to enforce loop charge neutrality) is imposed
to avoid numerical issues. Such an assumption is partly supported by the MD
simulations of loop formation in CeO2 aidhy2011comparison . However, the
diffusivity of oxygen interstitials is much higher than that of Th
interstitials colbourn1983calculated ; he2022dislocation  and the formation of
an initial loop may involve reconfiguration of small interstitial clusters
jiang2022unraveling . Hence, the dislocation loop can be off-stochiometric. To
quantify how strongly these oxygen atoms bind to the loop, FIG. 1c shows a
heatmap of the individual oxygen binding energy in a 10 nm diameter faulted
loop by calculating the energy difference of the loop with respect to placing
the loop constituent oxygen atom in the far field. The core oxygen atoms are
much more strongly bonded than the periphery oxygen atoms due to the local
lattice distortion at the dislocation core. As strong pipe diffusion was
identified in UO2 murphy2014pipe , it is possible that under irradiation
conditions, the loop tends to be oxygen-rich via trapping oxygen
interstitials. This may reduce the energy barrier for loop unfaulting due to
instability from charged loops, which is implied in the electron irradiation
experiments of ion-irradiated ZrO2, where fast loop growth and disappearance
were observed when only oxygen defects are created ryazanov2002growth ;
yasuda2003radiation .

The elastic strain field stukowski2012elastic  is computed near the loop as
shown in FIG. 2a-b for volumetric strain ($\epsilon_{V}$) and shear strain
($\epsilon_{YZ}$), respectively. As expected, the strain is exceptionally high
near the core region. The varying field of shear strain parallel to the plane
indicates a possible double-lobe contrast under TEM, as similarly indicated in
reference ryazanov2002growth  with diminished shear strain near the center. To
reveal the effect of strain field on atom kinetics, system annealing at 1500 K
and zero external pressure is performed. FIG. 2c demonstrates a heatmap of
oxygen displacements (Th atoms are removed due to negligible displacement)
after 400 ps. Large displacement of oxygen atoms mainly occurs at the loop
line, consistent with the previous understanding of pipe diffusion in ThO2 and
UO2 murch1987oxygen ; murphy2014pipe . Given the unstable state of oxygen
atoms at the dislocation core, it is expected that Shockley partials are
easier to nucleate at the loop circumference than within the loop.

Figure 2: Distribution of volumetric strain $\epsilon_{V}$ (a), shear strain
$\epsilon_{YZ}$ (b). (c) shows the oxygen displacement $|\Delta r|$ after 400
ps relaxation at 1500 K.

The unfaulting process can be in principle initiated due to shear stress or
thermal vibration kadoyoshi2007molecular ; chen2021comprehensive . Hence, two
methods have been attempted to initiate the transformation, i.e., shear and
high temperature. The shear is achieved by static shear of the simulation box,
via iterative i) rigid 0.1 Å displacement of the top layer along $y$-axis with
the bottom layer fixed, and ii) subsequent energy minimization. A sudden
unfaulting due to mechanical instability occurs when the simulation box is
sheared to 8.5 % as a result of lattice translation, without any indication of
Shockley partial nucleation (see Supplementary video 1). This phenomenon is
consistent with previous simulations of loop unfaulting in FCC Cu
kadoyoshi2007molecular . It suggests that thermal activation is necessary to
initiate the unfaulting process. Hence, the subsequent simulations are run at
1,500 K at zero external pressure to capture the unfaulting process. For a
reasonable computational cost, the initial configuration is slightly sheared
(2.68 %) along the $y$-axis, as this effectively reduces the energy barrier
for the nucleation of Shockley partial at a specific temperature. For
reference, high-temperature annealing up to 2,500 K of non-sheared structure
does not exhibit the unfaulting process on the nanosecond timescale. FIG. 3
demonstrates the complete unfaulting process within a span of 55 ps: Shockley
partial is nucleated from the loop line, and quickly sweeps across the fault
plane loop, which ultimately converts the original loop to a perfect one (see
Supplementary video 2). In another case with a 5.15 % sheared simulation cell,
a more complex unfaulting process is revealed, which involves the nucleation
of multiple Shockley partials along the loop line (FIG. 4). The propagation of
the arc segments contributes to the unfaulting simultaneously, leading to
reactions between Shockley partials within the loop plane. Since the nucleated
Shockley partials are not necessarily of the same type within the (111) plane,
the reaction can cause defect formation from local lattice incompatibility
(see Supplementary video 3). At the end of the unfaulting, the defects
gradually disappear due to atom diffusion. In both unfaulting patterns, by
analyzing the faults in both the Th and O layers, two common points are
identified: i) the extrinsic SF is removed in a single sweep of Shockley
partial; and ii) Both Th and O sublattice stacking faults are removed in a
synchronous manner.

Figure 3: Loop unfaulting process with single Shockley partial nucleation
site. Top row: dislocation types, middle row: Th stacking fault, bottom row: O
stacking fault. Figure 4: Loop unfaulting process with multiple Shockley
partial nucleation sites. Top row: dislocation types, middle row: Th stacking
fault, bottom row: O stacking fault.

We note some important observations of loop unfaulting behavior in this
system. First, it is feasible to have both single and multiple partial
nucleation. The large-sized loop utilized in the current work provides a
strong driving force to unfaulting. Such nucleation behavior is ultimately
determined by the local stress and temperature-activated atomic events. It has
also been reported in the unfaulting of both vacancy and interstitial faulted
loops in Cu, Al, and Ni kubota2005transition ; chen2021comprehensive ;
chen2022combined . Second, unlike the commonly proposed mechanisms for loop
unfaulting where a Shockley partial loop can be nucleated within the loop to
explain experimentally observed loop contrast chen2021comprehensive ;
chen2022combined , only Shockley partial nucleation at the loop perimeter is
observed in this system; this is presumably due to that the strong strain
field and oxygen mobility near the dislocation core (FIG. 1) can facilitate
the partial nucleation upon thermal fluctuation. Third, we have not observed
two-stage unfaulting as reported in metals kubota2005transition ;
kadoyoshi2007molecular ; chen2022combined , where a first set of Shockley
partials sweeps the extrinsic stacking fault, converting it into an intrinsic
stacking fault and a second set removes the remaining fault. To explain the
non-existence of any intrinsic stacking fault ribbon, it is only possible that
a joint Shockley partial pair sweep the fault simultaneously, with one
shearing the above and the other shearing the below atoms, respectively.
Hence, the total reaction is written as,

$\displaystyle 1/6[\bar{1}2\bar{1}]$
$\displaystyle+1/6[2\bar{1}\bar{1}]+1/3[111]\rightarrow 1/2[110]$
$\displaystyle A\delta$ $\displaystyle+B\delta+D\delta\rightarrow
DC~{}(\text{Thompson tetrahedron notation})$

Herein, $A\delta+B\delta\rightarrow\delta C(1/6[11\bar{2}])$, in the direction
of applied shear. The obtuse angle between the two partials suggests
attractive interaction, which can justify the close range of two joint
partials. Based on energetic estimation proportional to $\mathbf{b}^{2}$,
$|A\delta|^{2}+|B\delta|^{2}>|\delta C|^{2}$, but overall
$|A\delta|^{2}+|B\delta|^{2}+|D\delta|^{2}>|DC|^{2}$, leading to energy
reduction.

To explain the concurrent removal of faults in both Th and O layers given the
sweep of a joint Shockley partial pair, the details of layer stacking are
presented in the following. The layer stacking sequence along [111] is
depicted in FIG. 5a for the perfect fluorite structure, where cation, anion,
and octahedral site layers are interleaved to form an FCC stacking sequence
le2016empirical . We denote the sequence of cations, …abc…, the anions …ABC…,
and the octahedral sites, …$\alpha\beta\gamma$…, and the combined sequence is
…aB$\gamma$AbC$\alpha$BcA$\beta$C…. To construct the faulted loop, three
layers (O-Th-O) are inserted into the octahedral site layer, forming
…aB$\gamma$AbC/BaC/BcA$\beta$C… (FIG. 5b), where the italic font indicates the
fault layers. Hence, the interstitial loop means extrinsic stacking fault for
cation layers and intrinsic stacking fault for anion layers in the
corresponding sublattices. Using the Thompson tetrahedron notation, on $(111)$
plane, only the following Shockley partials can be accommodated, i.e.,
$A\delta$, $\delta A$, $B\delta$, $\delta B$, $C\delta$, and $\delta C$. To
align with simulation observations, with respect to the fault cation layer
$a$, layers above and below $a$ must be simultaneously sheared by a Shockley
partial, $1/6[\bar{1}2\bar{1}]$ $(A\delta)$ and $-1/6[2\bar{1}\bar{1}]$
$(-B\delta$, note the sign comes from the inverse view for the above and
below), respectively, resulting in …bC$\alpha$BcACaBAbC$\alpha$B… (FIG. 5c).
Therefore, with respect to the layer $a$, the relative displacement of the
adjacent layers (both cation and anion layers) should be around
$a_{0}/\sqrt{6}$, and diminishes as distance increases; this is confirmed by
the atom displacement field after unfaulting (FIG. 5d shows the displacement
along $y$-axis).

The O layers are generally commensurate with Th layers in displacement,
nevertheless, large displacement for O atoms close to the whole loop plane is
also observed (not shown in FIG. 5d). It is attributed to the pipe diffusion
mechanism as the Shockley partial sweeps across the plane. Similar to
dislocation nucleation, the high O diffusivity may reduce the energy barrier
of dislocation migration as it can induce local off-stoichiometry, causing a
metastable state of dislocation segments. After loop unfaulting, the perfect
loop shows high mobility along the Burgers’ vector direction (FIG. 5e) given
the applied shear strain. Under irradiation conditions, these perfect loops
can react with other extended defects to become dislocation networks.
Importantly, defective cation and anion atoms are found only exist at the
dislocation core of the perfect dislocation loops. Hence, unfaulting can be of
significance to controlling the total defect accumulation in ThO2. Similar to
FCC metals, a thermodynamic critical size of faulted loops should exist in
fluoride oxides, on the order of a few nanometers based on the formation
energy of loops he2022dislocation . Such critical size is generally smaller
than that of metals due to its high stacking fault energy (e.g., 1.76 J/m2 in
ThO2 he2022dislocation , 1.808 J/m2 in UO2 le2016empirical , and 78.5 mJ/m2 in
Cu heino1999stacking ). However, due to the nucleation barrier of Shockley
partials and interaction with other radiation-induced defects, this
thermodynamic critical size may not be particularly meaningful. For example,
in heavy ion irradiated CeO2, no perfect loop was observed for loops up to
20-30 nm in diameter, and in ThO2 and UO2 with various irradiation types,
faulted loops are frequently observed beyond a few nanometer sizes
he2022dislocation ; whapham1965radiation ; he2013situ ; onofri2017full .
Finally, it is worth noting that the unfaulting, specifically partial
nucleation, can be strongly affected by the radiation condition. For example,
only faulted loops are identified in proton-irradiated ThO2 bawane2021tem
while both types of loops are found in Kr-irradiated ThO2 he2022dislocation ;
this contrast could suggest that pronounced damage cascade disturbance assists
in the partial nucleation, which warrants further investigation.

Figure 5: Stacking sequences for perfect lattice (a), with faulted loop (b),
and with unfaulted loop (c). (d) shows a heatmap of displacement along
$y$-axis ($dY$) after loop unfaulting. (e) shows loop migration over $\sim$45
ps after loop unfaulting.

In summary, with large-size stoichiometric circular Frank loops, we have
identified two unfaulting patterns involving single and multiple Shockley
partial pairs under combined shear strain and high temperature to explain the
Frank loop unfaulting in fluorite ThO2. Each Shockley pair is jointly moving
in a synchronous manner to transform the stacking faults into perfect
stacking. Cation and anion layers are commensurate during loop unfaulting and
perfect loop migration. The high mobility of anions along dislocation line
likely reduces the activation barriers of dislocation activities, which
remains to be confirmed.

## Acknowledgments

This work is supported by the Center for Thermal Energy Transport under
Irradiation, an Energy Frontier Research Center funded by the U.S. Department
of Energy, Office of Science, United States, Office of Basic Energy Sciences,
and The Pennsylvania State University.

## References

  * [1] Dongyue Chen, Kenta Murakami, Kenji Dohi, Kenji Nishida, Zhengcao Li, and Naoto Sekimura.  The effects of loop size on the unfaulting of frank loops in heavy ion irradiation.  Journal of Nuclear Materials, 529:151942, 2020.
  * [2] Alison Kubota and WG Wolfer.  Transition pathways in the unfaulting of dislocation loops.  Materials Science and Engineering: A, 400:362–365, 2005.
  * [3] Lingfeng He, Tiankai Yao, Kaustubh Bawane, Miaomiao Jin, Chao Jiang, Xiang Liu, Wei-Ying Chen, J Matthew Mann, David H Hurley, Jian Gan, et al.  Dislocation loop evolution in kr-irradiated tho2.  Journal of the American Ceramic Society, 105(8):5419–5435, 2022\.
  * [4] Wei-Ying Chen, Jianguo Wen, Marquis A Kirk, Yinbin Miao, Bei Ye, Brian R Kleinfeldt, Aaron J Oaks, and James F Stubbins.  Characterization of dislocation loops in ceo2 irradiated with high energy krypton and xenon.  Philosophical Magazine, 93(36):4569–4581, 2013.
  * [5] Bernd Baufeld, Dietmar Baither, Ulrich Messerschmidt, Martin Bartsch, and Ingolf Merkel.  In situ study on the generation of radiation damage in cubic-zirconia in the high-voltage electron microscope.  Journal of the American Ceramic Society, 76(12):3163–3166, 1993\.
  * [6] K Yasunaga, Kazuhiro Yasuda, Syo Matsumura, and T Sonoda.  Electron energy-dependent formation of dislocation loops in ceo2.  Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 266(12-13):2877–2881, 2008.
  * [7] AI Ryazanov, Kazuhiro Yasuda, C Kinoshita, and AV Klaptsov.  Growth and instability of charged dislocation loops under irradiation in ceramic materials.  Journal of nuclear materials, 307:918–923, 2002.
  * [8] Derek Hull and David J Bacon.  Introduction to dislocations.  Butterworth-Heinemann, 2001.
  * [9] Tomoko Kadoyoshi, Hideo Kaburaki, Futoshi Shimizu, Hajime Kimizuka, Shiro Jitsukawa, and Ju Li.  Molecular dynamics study on the formation of stacking fault tetrahedra and unfaulting of frank loops in fcc metals.  Acta materialia, 55(9):3073–3080, 2007.
  * [10] Cheng Chen, Jing Zhang, and Jun Song.  Comprehensive study of vacancy frank loop unfaulting: atomistic simulations and predictive model.  Acta Materialia, 208:116745, 2021.
  * [11] C Dai, Q Wang, P Saidi, B Langelier, CD Judge, MR Daymond, and MA Mattucci.  Atomistic structure and thermal stability of dislocation loops, stacking fault tetrahedra, and voids in face-centered cubic fe.  Journal of Nuclear Materials, 563:153636, 2022.
  * [12] Cheng Chen and Jun Song.  A combined atomistic-continuum study on the unfaulting of single and multi-layer interstitial dislocation loops in irradiated fcc and hcp metals.  International Journal of Plasticity, 152:103231, 2022.
  * [13] Steve Plimpton, Paul Crozier, and Aidan Thompson.  LAMMPS-large-scale atomic/molecular massively parallel simulator.  Journal of Computational Physics, 18, 2007.
  * [14] MWD Cooper, MJD Rushton, and RW Grimes.  A many-body potential approach to modelling the thermomechanical properties of actinide oxides.  Journal of Physics: Condensed Matter, 26(10):105401, 2014.
  * [15] Shuichi Nosé.  A unified formulation of the constant temperature molecular dynamics methods.  The Journal of Chemical Physics, 81(1):511–519, 1984.
  * [16] William G Hoover.  Canonical dynamics: Equilibrium phase-space distributions.  Physical Review A, 31(3):1695, 1985.
  * [17] A Chartier, Claire Onofri, L Van Brutzel, Ch Sabathier, O Dorosh, and J Jagielski.  Early stages of irradiation induced dislocations in urania.  Applied Physics Letters, 109(18):181902, 2016.
  * [18] Arno Le Prioux, Paul Fossati, Serge Maillard, Thomas Jourdan, and Philippe Maugis.  Empirical potential simulations of interstitial dislocation loops in uranium dioxide.  Journal of Nuclear Materials, 479:576–584, 2016.
  * [19] Dilpuneet S Aidhy, Dieter Wolf, and Anter El-Azab.  Comparison of point-defect clustering in irradiated ceo2 and uo2: A unified view from molecular dynamics simulations and experiments.  Scripta Materialia, 65(10):867–870, 2011.
  * [20] Ling-Feng He, Mahima Gupta, Clarissa A Yablinsky, Jian Gan, Marquis A Kirk, Xian-Ming Bai, Janne Pakarinen, and Todd R Allen.  In situ tem observation of dislocation evolution in kr-irradiated uo2 single crystal.  Journal of Nuclear Materials, 443(1-3):71–77, 2013.
  * [21] Alexander Stukowski.  Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool.  Modelling and Simulation in Materials Science and Engineering, 18(1):015012, 2009.
  * [22] EA Colbourn and WC Mackrodt.  The calculated defect structure of thoria.  Journal of nuclear materials, 118(1):50–59, 1983.
  * [23] Chao Jiang, Lingfeng He, Cody A Dennett, Marat Khafizov, J Matthew Mann, and David H Hurley.  Unraveling small-scale defects in irradiated tho2 using kinetic monte carlo simulations.  Scripta Materialia, 214:114684, 2022.
  * [24] Samuel T Murphy, Eleanor E Jay, and Robin W Grimes.  Pipe diffusion at dislocations in uo2.  Journal of Nuclear Materials, 447(1-3):143–149, 2014.
  * [25] K Yasuda, C Kinoshita, S Matsumura, and AI Ryazanov.  Radiation-induced defect clusters in fully stabilized zirconia irradiated with ions and/or electrons.  Journal of Nuclear Materials, 319:74–80, 2003.
  * [26] Alexander Stukowski and A Arsenlis.  On the elastic–plastic decomposition of crystal deformation at the atomic scale.  Modelling and Simulation in Materials Science and Engineering, 20(3):035012, 2012.
  * [27] GE Murch and C Richard A Catlow.  Oxygen diffusion in uo 2, tho 2 and puo 2. a review.  Journal of the Chemical Society, Faraday Transactions 2: Molecular and Chemical Physics, 83(7):1157–1169, 1987.
  * [28] Pekka Heino, L Perondi, Kimmo Kaski, and E Ristolainen.  Stacking-fault energy of copper from molecular-dynamics simulations.  Physical Review B, 60(21):14625, 1999.
  * [29] AD Whapham and BE Sheldon.  Radiation damage in uranium dioxide.  Philosophical Magazine, 12(120):1179–1192, 1965.
  * [30] C Onofri, M Legros, J Léchelle, H Palancher, C Baumier, C Bachelet, and C Sabathier.  Full characterization of dislocations in ion-irradiated polycrystalline uo2.  Journal of Nuclear Materials, 494:252–259, 2017.
  * [31] Kaustubh Bawane, Xiang Liu, Tiankai Yao, Marat Khafizov, Aaron French, J Matthew Mann, Lin Shao, Jian Gan, David H Hurley, and Lingfeng He.  Tem characterization of dislocation loops in proton irradiated single crystal tho2.  Journal of Nuclear Materials, 552:152998, 2021.


Current address: ]Department of Condensed Matter Physics, Weizmann Institute
of Science, Rehovot, Israel

# Magnetic interactions of 4$f$ electrons in the topological insulator
chalcogenide Bi2Se3

J. C. Souza<EMAIL_ADDRESS>[  Instituto de Física “Gleb Wataghin”,
UNICAMP, 13083-859, Campinas, SP, Brazil  M. Carlone  POSMAT, Faculdade de
Ciências, UNESP, C.P. 473, 17033-360, Bauru, SP, Brazil  G. G. Lesseux  1\.
Physikalisches Institut, Universität Stuttgart, D-70569, Germany  H. B. Pizzi
Instituto de Física “Gleb Wataghin”, UNICAMP, 13083-859, Campinas, SP, Brazil
G. S. Freitas  Instituto de Física “Gleb Wataghin”, UNICAMP, 13083-859,
Campinas, SP, Brazil  Los Alamos National Laboratory, Los Alamos, New Mexico
87545, USA  R. R. Urbano  Instituto de Física “Gleb Wataghin”, UNICAMP,
13083-859, Campinas, SP, Brazil  P. A. Venegas  Departamento de Física,
Faculdade de Ciências, UNESP, C.P. 473, 17033-360, Bauru, SP, Brazil  P. G.
Pagliuso  Instituto de Física “Gleb Wataghin”, UNICAMP, 13083-859, Campinas,
SP, Brazil  Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA

###### Abstract

The gap opening mechanism of a topological insulator, the quantum anomalous
Hall effect and the axion physics are still pressing open questions and a
microscopic viewpoint to further understand the role of magnetism in topology
is highly desirable. In this work we have performed a microscopic
investigation, by means of electron spin resonance (ESR) along with
complementary bulk measurements, on the chalcogenide (Bi1-xGdx)2Se3 ($x$ = 0,
0.001, 0.002 and 0.006). Our analysis of the Gd3+ spin dynamics reveal no
significant change of the Fermi surface as a function of Gd3+ concentration,
which indicates that the 4$f$ magnetism is different from the non-local
effects induced by transition metals ($d$ electrons) substitutions.
Additionally, we observe an unusual evolution of the Gd3+ ESR spectra as a
function of the applied magnetic field, which we discuss considering the
magnetic interaction between Gd3+ 4$f$ electrons and impurity centers such as
Se vacancies. This interaction would give rise to a local weak
antilocalization effect surrounding the Gd3+ ions. Such mechanism is
observable due to particular details of the Gd3+ 4$f$ electrons magnetism in
this system compared to $d$ electrons. Our work points out that rare earth
substitutions in this model topological insulator is a promising path to
explore the axion insulating systems.

###### pacs:

76.30.-v, 71.20.Lp

††preprint: APS/123-QED

##  I I. Introduction

The application of topology in condensed matter physics has been responsible
to unveil new quantum phases of matter, the topological insulators (TIs) in
two [1, 2, 3] and three dimensions [4, 3] being the most prominent and heavily
explored systems. The insulating bulk with protected spin-polarized gapless
surface states is a highly attractive characteristic of TIs [4, 5]. Another
interesting property rises from the interplay between topology and magnetism,
where the time reversal symmetry breaking could result in axion physics and
the quantum anomalous hall effect (QAHE) [6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17].

The (Bi,Sb)2(Te,Se)3 chalcogenides have been established as the model systems
of three dimensional TIs due to the single Dirac cone near the Fermi level
[18, 19, 20]. Naturally, such systems have been able to host the interplay
between magnetism and topology [7, 8, 9]. The first realization of the QAHE
was achieved in (V,Cr):(Bi,Sb)2Te3 thin films, however it is possible to
obtain a fully quantized Hall conductivity only at milikelvin temperatures [7,
8, 9, 10]. This limitation may be due to the presence of thermally activated
bulk carriers, which cause a breakdown of the QAHE [21]. One remaining
question relies on the role of the underlying magnetic coupling mechanism.
More specifically, it is imperative to understand how the QAHE can be affected
by more complex magnetic interactions which goes beyond the so-called van
Vleck mechanism. This question has recently been explored by electron spin
resonance (ESR), x-ray absorption and resonant photoelectron spectroscopy [22,
23, 24, 25, 26, 27, 28, 29]. Indeed the $p-d$ hybridization, along with the
$d$ states occupation and the consequent Sb and Te polarization, near the
Fermi level seems to play an important role to the magnetism of substituted
chalcogenides, which goes beyond the van Vleck mechanism [25, 30, 26, 29].
Exploring heterostructures is also a promising, yet challenging, path to
achieve higher temperatures in the QAHE [17, 31, 32, 11, 33, 34, 35].

One little explored but promising route to understand the role of the
magnetism and its influence in the bulk is to investigate the magnetism of
4$f$ electrons. For instance, Sm-substituted chalcogenides could display
higher order topological insulators phases with chiral hinge states [36, 12].
Europium substituted chalcogenides shows antiferromagnetic correlations,
however Eu has a 2+ valence, while most of the rare-earths shows a 3+ state
[37]. Angle resolved photoemission spectroscopy studies show the robustness of
the surface states to a Gd3+ concentration of $\sim$ 0.1 in Bi2Te3 and TlBiSe2
[38, 39]. Therefore, further exploring 4$f$-substituted systems can be highly
interesting to understand the gap opening mechanism and the axion insulating
phase in these model systems [40, 41, 42, 39, 43].

In this context, Gd3+, which also induces antiferromagnetic correlations [44,
45, 42], is an ideal testbed due to the stable valence and the weak
crystalline electrical field (CEF) effects. Although Gd3+-substituted
chalcogenides have been explored in macroscopic [44, 45] and previous ESR
studies [44, 46, 47, 48], a detailed investigation of the 4$f$ local magnetic
effects and the spin dynamics induced by these substitutions, as well as a
comparison to more traditional $d$ systems is missing [27, 23].

In this work we locally explore Gd3+-substituted Bi2Se3 using ESR at different
frequencies ($\nu$ = 9.5 and 34 GHz). We show that the introduction of Gd3+
ions does not alter the carriers near the Fermi surface, which are mainly $p$
states. This is a different mechanism from the substitution using ions where
the magnetism comes from $d$ states. Such difference is also manifested into
the Gd3+ spin dynamics and in the macroscopic properties of the system.
Additionally, the Gd3+ ESR data show an unusual evolution as a function of the
applied magnetic field. While for lower field we obtain a Gd3+ response with
resolved fine structure, which is more consistent with Gd3+ ions in an
insulating environment, at higher fields the system tends to behave as a
conventional metal, revealing a single additional line with collapsed fine
structure. We discuss the evolution of the Gd3+ local environment under the
light of a possible local weak antilocalization (WAL) effect [49, 50, 51, 52],
which is a product of the interplay between strong spin-orbit coupling and the
interaction between 4$f$ local moments and impurity centers such as Se
vacancies. Our results shed light on the microscopic mechanism involving the
introduction of 4$f$ electrons into the chalcogenides.

##  II II. Methods

The chalcogenides have a rhombohedral crystal structure. Single crystalline
samples of (Bi1-xGdx)2Se3 ($x$ = 0, 0.001, 0.002 and 0.006) were grown by the
stoichiometric melting technique. High purity Bi, Gd and Se elements were put
inside of an alumina crucible with the ratio of (2 - $2x$):$2x$:3. The
crucible was vacuum sealed in a quartz tube, heated to 800∘C for 72 hours and
cooled down to room temperature at 2∘C per hour. We cut the crystals in a
rectangular-like shape. Typical sample sizes were 0.5 x 2 x 0.3 mm3. The
structure and phase purity were checked by x-ray powder diffraction using a
commercial diffractometer (Cu-K$\alpha$), where we confirmed the single phase
nature of our samples. We have performed elemental analysis using energy and
wavelenght dispersive spectroscopy, which revealed a small, but measurable,
amount of Se vacancies ($\sim$ 0.1). Magnetic susceptibility measurements were
performed in a commercial SQUID magnetometer. Specific-heat measurements were
done in a commercial small-mass calorimeter system. Electrical resistivity
data were acquired using a four-probe configuration with a $dc$ resistance
bridge. ESR measurements were performed on single crystals in X- and Q-band
($\nu$ = 9.5 and 34 GHz respectively) spectrometer equipped with a goniometer
and a He-flow cryostat in the temperature range of 4 K $\leq$ $T$ $\leq$ 300 K
at low power $P$ $\leq$ 2 mW. The ESR spectra were analyzed using the software
SPEKTROLYST.

##  III III. Results

Figure 1: (Top panel) Magnetic susceptibility for $H$ applied parallel to the
[001] direction, (middle panel) longitudinal resistivity and (lower panel)
specific heat as a function of temperature for (Bi1-xGdx)2Se3. The insets show
the (top panel) susceptibility at high temperatures, (middle panel) Hall
resistivity as a function of applied magnetic field for different temperatures
for $x$ = 0.002, and (bottom panel) the specific heat divided by temperature
as a function of $T^{2}$ at low temperatures.

Figure 1 summarizes the macroscopic properties of (Bi1-xGdx)2Se3. The top
panel shows the magnetic susceptibility as a function of temperature. From a
Curie-Weiss fitting $\chi=\chi_{0}+C/(T-\theta$), where $\chi_{0}$ is the
$T$-independent term, $C$ is a constant and $\theta$ the Curie-Weiss
temperature, at high temperatures (150 K $\leq$ $T$ $\leq$ 300 K) we can
estimate the concentration of Gd3+ in (Bi1-xGdx)2Se3. Assuming 7.94
$\mu_{B}$/Gd, we obtained $x$ $\approx$ 0.001, 0.002 and 0.006 with a $\theta$
= -1(5) K for all samples. Taking into account the nuclear diamagnetism, we
can also estimate the Pauli magnetic susceptibility $\chi_{p}$ = 20(10)
$\mu$emu/mol Oe for all samples.

The middle panel shows the longitudinal resistivity $\rho_{xx}$ as a function
of temperature for (Bi1-xGdx)2Se3. We obtain a metallic-like behavior, which
is expected due to the presence of Se vacancies in single crystals [50, 53,
54, 55]. Nonetheless, we only observe small and not systematic differences
between samples, which can be attributed to a small variation of Se vacancies
from crystal to crystal. As a result, the residual resistivity $\rho_{0}$,
which can be associated with disorder in the system, can differ between
samples from different batches. The lack of a systematic change indicates that
Gd3+ substitution at the Bi site does not introduce carriers into the system,
and presumably its 4$f$ electrons have no relevant role in the Fermi surface.
This is in contrast to substitutions using transition metal ions where the
magnetism originates from $d$ orbitals [56, 23, 29, 57]. The absence of
Gd3+-introduced carriers into the system as a function of Gd3+ concentration
is also supported by the Hall resistivity ($\rho_{xy}$), which is Gd3+
concentration and temperature independents, as shown in the inset of the
middle panel. The Hall response to the applied magnetic field is linear and
positive, consistent with the transport properties being dominated by a single
band of holes ($p$-type). From a linear fit ($\rho_{xy}$ = 1/ne, where $e$ is
the electron charge) we can extract a carrier density of $n_{h}$ = 5(3).1018
h/cm3.

The bottom panel shows the specific heat divided by the temperature
$c_{p}$/$T$ as a function of temperature for (Bi1-xGdx)2Se3. There is a slight
difference between different Gd3+ concentrations, which reinforces the lack of
change in the Fermi level due to the Gd3+ substitution. Performing linear fits
for $c_{p}$/$T$ as a function of $T^{2}$, top inset of the bottom panel, it is
possible to extract the Sommerfeld coefficient $\gamma$ = 0.8(3) mJ/mol K2 for
all samples. For a free conduction electron gas model $\gamma=(2/3)\pi
k_{B}^{2}\eta(E_{F})$, where $k_{B}^{2}$ is the Boltzmann constant and
$\eta(E_{F})$ the density of states (DOS) at the Fermi level per spin. We
extract $\eta(E_{F})$ = 0.16(6) states/eV mol spin for all samples. We can
further analyze the role of electron-electron ($ee$) interactions in the bulk
by comparing the estimated Pauli susceptibility ($\chi_{p}^{theor}$ =
2$\mu_{B}^{2}\eta(E_{F})$, where $\mu_{B}$ is the Bohr magneton) with the
experimental value. We obtain $\chi_{p}^{theor}$ = 15 $\mu$emu/mol Oe. These
results indicate that the DOS at the Fermi level is not affected as a function
of Gd3+ concentration and, if any, the role of $ee$ interactions is negligible
in the bulk. It is noteworthy that recent results point out that $ee$
correlations do play a role in the surface states of chalcogenides [58, 59,
60].

Figure 2: (a) X-band ($\nu$ = 9.5 GHz) Gd3+ ESR spectra for
(Bi0.994Gd0.006)2Se3 and (b) its resonance field angle dependence at $T$ = 4 K
with the applied magnetic field $H$ parallel to the [001] direction. In (b)
the field is rotated towards the $ab$ plane. (c) Q-band ($\nu$ = 34 Ghz) Gd3+
ESR spectra. The blue solid lines are simulations considering the relative
intensities for each of the seven Gd3+ fine structure lines. The magenta solid
lines are fits considering the hexagonal CEF Hamiltonian described into the
text.

As mentioned earlier, Gd3+ is a $S$-ion ($L$ = 0, $S$ = 7/2). As such, the CEF
effects appears only due to a small mixing of excited states in the ground
state, resulting in a intermediate coupling [61, 62]. Due to this small
correction, the CEF splitting is of the order of the Zeeman energy and we are
able to observe the Gd3+ fine structure, with a selection rule of $\Delta
m_{s}$ = $\pm$ 1 [61, 62]. Another consequence of the $S$-ion nature is that
the $g$-value of an isolated Gd3+ ion is independent of the symmetry of the
matrix, which makes the analysis of the Gd3+ spin dynamics more
straightforward [61, 62]. With that in mind, in order to have a microscopic
insight about the local effects of Gd3+ ions, we have performed ESR in two
different frequencies. Figure 2 shows the Gd3+ ESR spectra at $T$ = 4 K for
$x$ = 0.006 with the applied field $H$ parallel to the [001] direction and its
angle dependence as an example. Fig. 2(a) shows the Gd3+ X-band ESR spectrum,
where we can observe the resolved Gd3+ fine structure displaying seven
transition lines [61, 62, 63] \- with their expected relative intensities.
Such observation excludes the possibility of Gd3+ ions having distinct sites
and/or being interstitial in Bi2Se3. It also corroborates that we only have
one phase in our crystals. Each individual resonance line can be described as
a Dysonian line shape, which is typical when the skin depth $\delta$ is
smaller than the sample size $d$ [64, 65]. Indeed, from the resistivity
measurements shown in Fig. 1(b), we can estimated for X-band $\delta$ $\sim$
11 $\mu$m ($\ll$ $d$ $\sim$ 300 $\mu$m for our samples) at $T$ = 4 K,
consistent with the Dysonian line shape. Each individual line, represented by
the power absorption derivative d$P$/d$H$ as a function of $H$, can be
described as an admixture of absorption and dispersion derivatives

$\frac{dP}{dH}\propto(1-\lambda)\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{1}{1+x^{2}}\right)+\lambda\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{x}{1+x^{2}}\right),$
(1)

where $\lambda$ is the asymmetric parameter of the line shape and $x$ =
$2(H-H_{r})/\Delta H$ [64], wherein $H_{r}$ is the Gd3+ resonance field and
$\Delta H$ the linewidth.

As already mentioned, the position of the Gd3+ fine structure is dependent on
the local symmetry of Gd3+ ions. In the case of the Bi2Se3, the local symmetry
is hexagonal, where the spin Hamiltonian is given by

$\begin{split}\mathcal{H}=&\frac{1}{3}b_{2}^{0}O_{2}^{0}+\frac{1}{60}\left(b_{4}^{0}O_{4}^{0}+b_{4}^{3}O_{4}^{3}\right)+\\\
&\frac{1}{1260}\left(b_{6}^{0}O_{6}^{0}+b_{6}^{3}O_{6}^{3}+b_{6}^{6}O_{6}^{6}\right)+g\mu_{B}\mathbf{H}\mathbf{S},\end{split}$
(2)

where $b_{n}^{m}$ are the $n$ order CEF parameters and $O_{m}^{n}$ the Stevens
operator. The Gd3+ fine structure resonance fields are going to be highly
angular dependent [61, 62]. Therefore, we fitted the angular dependence of the
Gd3+ fine structure resonance fields with Eq. 2, as shown by the solid magenta
lines in Fig. 2(b). The best parameters obtained were $b_{2}^{0}$ = 37.5(5)
Oe, $b_{4}^{3}$ = 0.03(2) Oe and $b_{6}^{6}$ = - 0.3(2) Oe. The values of
$b_{4}^{0}$, $b_{6}^{0}$ and $b_{6}^{3}$ were negligible. The dominant term,
which is the axial $b_{2}^{0}$, is consistent with previous CEF studies of
Gd3+ in Bi2Se3 [47]. Albeit the contribution from smaller terms may change
from this previous study to our result, it is important to observe that the
absolute signs of $b_{2}$ $>$ 0, $b_{4}$ $>$ 0 and $b_{6}$ $<$ 0 are in
agreement between both studies [47]. The change between which terms are
negligible from one to another may be due to those terms being much smaller
than the dominant one, and we can obtain a local minimum different for each
rotation. Nonetheless, our obtained CEF parameters are in accordance with
[47]. The blue solid line is a simulation considering the value of $b_{2}^{0}$
of this hexagonal CEF Hamiltonian.

At this point it is interesting to compare our X-band Gd3+ spectrum in Bi2Se3
with previous reports of Gd3+ and Mn2+ in Bi2Te3 [23, 46, 44]. While both (X-
and Q-bands) Gd3+ spectra show resolved fine structure, the same does not
happen for Mn2+, in which only one resonance is observed even for the smallest
Mn2+ concentration ($\sim$ 0.005) [23]. Mn2+ is a $S$ = 5/2 ion, therefore
there is an exchange narrowing of the fine structure going from $d$ to $f$
states, which reflects different substitution effects. The exchange narrowing
will result from the large exchange interaction from $d$ states with the $p$
polarized states in these systems [56, 23, 29, 66, 67].

Fig. 2(c) shows the Gd3+ Q-band ESR spectrum at $T$ = 4 K. The blue solid line
is a simulation of the same CEF scheme of the X-band data, normalized by the
intensity of the Gd3+ fine structure. There is a clear difference when
compared to the X-band data, with a stronger - 1/2 $\leftrightarrow$ 1/2
transition, which we will further call Gd3+ central line, at the top of the
seven fine structure lines. Such coexistence can be linked to two different
environments for Gd3+ ions distributed along the sample. Although the
distribution of the Gd3+ ions appears to be structurally homogeneous, due to
the X-band data, the electronic environment around the localized moments seems
to show an interesting evolution as a function of the applied magnetic field.
Such effect may be connected with the Se vacancies distribution, as we discuss
below analyzing the Gd3+ spin dynamics in more details.

Figure 3: (a) Gd3+ spectrum splitting and the axial term $b_{2}^{0}$, (b)
$\Delta H$ and (c) $g$-shift as a function of temperature for
(Bi0.994Gd0.006)2Se3 for the applied magnetic field parallel to the [001]
direction. $\Delta H$ and $g$-shift were extracted from the ESR central line.
The red dashed lines in (b) represent the linear fits used to extract the
Korringa rates. We used $g_{theor}$ = 1.993(1) in order to calculate the
$g$-shift for all samples.

The Gd3+ ESR spectrum temperature evolution and the spin dynamics can bring up
valuable information about the introduction of 4$f$ states in Bi2Se3 and their
field dependence. Figure 3 summarizes this temperature evolution. The ESR fine
structure split, defined by the difference between the resonance field $H_{r}$
of the $\mp$ 7/2 $\leftrightarrow$ $\mp$ 5/2 transitions ($H_{r7}$ \-
$H_{r1}$), is shown in Fig. 3(a). Since the axial term of the CEF parameters
is the most dominant one, $b_{2}^{0}$ is going to have to be directly
proportional to the CEF splitting. We also show the value of $b_{2}^{0}$ for
each CEF splitting in Fig. 3(a). The Gd3+ splitting is nearly constant
($\Delta_{H_{r7}-H_{r1}}^{4-50K}$/$T$ $\leq$ 0.4 Oe/K) until $\sim$ 50 K for
both bands. At higher temperatures we observe a more significant
($\Delta_{H_{r7}-H_{r1}}^{50-300K}$/$T$ $\sim$ 1 Oe/K) and systematic
reduction of the fine structure splitting. This narrowing can be related with
the interaction between carriers or, more likely, due to changes in the CEF
effects related to the thermal expansion of the compound [47]. Indeed,
previous neutron diffraction and pair density function analysis have shown a
local anharmonic thermal expansion in Bi2Se3 [68]. In particular, the two Se
sites have a distinct thermal expansion in reference to the Bi site [68].
Therefore, this anharmonic expansion should affect the CEF effects in the Bi
site and, consequently, should influence the Gd3+ CEF splitting. Our results
are in agreement with previous results [47, 48] and shows the similarity of
the fine structure between X- and Q-bands spectra.

From the analysis using an admixture of absorption and dispersion derivatives
using Eq. 1 one can extract the Gd3+ $\Delta H$ and $H_{rc}$ of the central
line. As a result, we can also obtain the experimental Gd3+ $g$-value
$g_{exp}$ = $h\nu/\mu_{B}H_{rc}$, where $h$ is the Planck constant and
$\mu_{B}$ the Bohr magneton. Figures 3(b) and 3(c) show the Gd3+ linewidth
$\Delta H$ of the central line and $g$-shift $\Delta g$ = $g_{exp}$ \-
$g_{theor}$ as a function of temperature. Here $g_{theor}$ = 1.993(1), which
is the Gd3+ $g$-value in an insulating matrix. The Gd3+ central line is a
reliable resonance to analyze the trend of the Gd3+ spin dynamics due to a
smaller influence of CEF effects in the $g$-value [62].

Regarding the $\Delta H$ of the Gd3+ central line, all the obtained data show
a linear increase at high temperatures ($T$ $\geq$ 80 K). We focused our
analysis on the high-$T$ region to avoid the contribution of possible
Gd3+-Gd3+ interactions at low-$T$, especially for samples with higher
concentrations and higher fields. Such linear increase can be attributed to a
relaxation process through the exchange interaction between the Gd3+ 4$f$
local moments and the carriers, which eventually results in the spin-flip
scattering of the latter. This spin-spin relaxation mechanism is known as
Korringa relaxation [61, 62, 63]. From a linear fit $\Delta H$ = $\Delta
H_{0}$ \+ $bT$ we can extract the Gd3+ residual linewidth $\Delta H_{0}$ and
the Korringa rate $b$. The results are summarized in Table 1. While $b$ is
related to the spin-spin relaxation, $\Delta H_{0}$ can be associated to,
e.g., disorder and sample homogeneity [61, 62, 63].

Table 1: ESR parameters extracted from the Gd3+ spin dynamics analysis for (Bi1-xGdx)2Se3. | X-band | Q-band
---|---|---
|  $\Delta H_{0}$ |  $b$ |  $\langle J_{fp}^{2}(\textbf{q})\rangle^{1/2}$ |  $J_{fp}(0)$ |  $\Delta H_{0}$ |  $b$ |  $\langle J_{fp}^{2}(\textbf{q})\rangle^{1/2}$ |  $J_{fp}(0)$
|  Oe |  Oe/K |  meV |  meV |  Oe |  Oe/K |  meV |  meV
$x$ = 0.006 |  37(2) |  0.015(5) |  $\sim$ 5 |  $\sim$ 80 |  38(2) |  0.020(5) |  $\sim$ 5 |  $\sim$ 20
$x$ = 0.002 |  29(2) |  0.016(5) |  $\sim$ 5 |  $\sim$ 80 |  36(2) |  0.017(5) |  $\sim$ 5 |  $\sim$ 20
$x$ = 0.001 |  28(2) |  0.013(5) |  $\sim$ 5 |  $\sim$ 80 |  31(2) |  0.015(5) |  $\sim$ 5 |  $\sim$ 20

Turning our attention to the Gd3+ X-band results, $b$ is, within experimental
uncertainty, the same for all concentrations, which indicates the absence of
exchange bottleneck effects [61, 62, 63]. $\Delta H_{0}$ increases
systemically as a function of Gd3+ concentration, which is expected due to the
increase of disorder in the system. At low temperatures we see a deviation of
the Korringa rate, which might be associated with the interaction between
local moments and impurity centers. As such, as mentioned earlier, a proper
analysis is to focus on the high temperature Gd3+ spin dynamics data. Looking
now to the field dependency, we still obtain the same Korringa rate for all
concentrations [Fig. 3(b)]; however $\Delta H_{0}$ increases systematically
when comparing X- and Q-bands results for the same concentration. This is an
indication of inhomogeneous broadening, which raises most likely from slightly
different CEF states around the Gd3+ sites [61, 62, 63]. Bi2Se3 is well known
to host Se vacancies intrinsically, which can show a small inhomogeneity
across the sample [50, 69, 70]. Moreover, previous nuclear magnetic resonance
measurements have shown that, in polycrystalline Bi2Se3, defect regions
segregate into domains [71].

Before analyzing the spin relaxation of our Gd3+ probe in more details, it is
instructive to first look to the $g$-shift as a function of temperature,
applied magnetic field and Gd3+ concentration - Figure 3(c). On one hand, $s$
and/or $d$ carriers have a ferromagnetic (atomic) interaction with 4$f$ local
moments, which produces a positive $g$-shift. On the other hand, $p$ and/or
$f$ carriers magnetic interaction with 4$f$ local moments occurs through the
so-called virtual bound states [72], which results in an antiferromagnetic
(covalent) interaction between them. In this second case, the result is a
negative $g$-shift. Therefore, the sign of the $g$-shift is crucial to have
information about the nature of the wave function of the carriers of the
system.

As such, an important result reported here is the negative $g$-shift,
confirming $p$ states as the main carriers near the Fermi level [18]. This
result is consistent with previous angle resolved spectroscopy studies [73,
74]. Looking more specifically to the X-band data, at low temperatures we can
see a systematic decrease of the $g$-shift for all concentrations. This is due
to an antiferromagnetic interaction, which is in accordance with the observed
AFM ground state in Eu2+ and Gd3+-substituted chalcogenides [75, 37]. However,
such decrease is observed even for samples with a concentration as low as $x$
= 0.001, further suggesting that such contribution comes from a possible
interaction between 4$f$ local moments and impurity centers, such as Se
vacancies, which also have a signature in the Gd3+ $\Delta H$ $T$-dependence -
Fig. 3 (b). Similar contributions recently observed in organic salts were also
interpreted to have origin in impurity centers [76, 77, 78, 79, 80]. At high
temperatures, the Gd3+ $\Delta g$ values are $T$-independent within
experimental uncertainty, which shows that dynamic effects are negligible into
the systems [81]. While the Gd3+ $g$-values for $x$ = 0.001 and 0.006 are
virtually identical, as expected, for $x$ = 0.002 we observe a subtle, but
systematic, increase of the $g$-value. Such small difference might also be
related with samples with $x$ = 0.002 showing a higher quantity of impurity
centers, such as Se vacancies. This interpretation is corroborated by a few
points: The first one is that the increase of $\Delta H_{0}$ for $x$ = 0.002
going from X- to Q-band is relatively larger when compared to other Gd3+
concentrations. Secondly, the deviation at low temperatures of the linear
increase in the Gd3+ $\Delta H$ are more pronounced for $x$ = 0.002 when
compared to other concentrations. Finally, the $g$-shift for all three
concentrations are the same for Q-band measurements, which indicates that
there is no change in $\eta(E_{F})$ as a function of Gd3+ concentration,
consistent with specific heat measurements.

Interestingly, there is also a decrease of the negative Gd3+ $g$-shift as a
function of the applied magnetic field for all concentrations of
(Bi1-xGdx)2Se3. In the absence of dynamic, bottleneck and multiple band
effects, $b$ and $\Delta g$ can be described as

$b=\frac{\pi k_{B}}{g\mu_{B}}\langle
J_{fp}^{2}(\textbf{q})\rangle\eta^{2}(E_{F})\frac{K(\alpha)}{(1-\alpha)^{2}},$
(3)

$\Delta g=J_{fp}(0)\,\frac{\eta(E_{F})}{(1-\alpha)},$ (4)

where $k_{B}$ is the Boltzmann constant, $J_{fp}(0)$ is the effective exchange
interaction between the Gd3+ local moments and the carriers for the momentum
transfer $q$ = 0, and $\langle J_{fp}^{2}(\textbf{q})\rangle^{1/2}$ is the
average of the exchange interaction with momentum transfer $q$ at the residual
Fermi surface [61, 62, 63]. Any possible relevant $ee$ correlations are taken
into account in the Stoner enhancement factor $(1-\alpha)^{-1}$ [82, 83] and
in the Korringa exchange factor $K(\alpha)$ [84, 85]. As already shown, the
$ee$ correlations in the bulk does not seem to play an important role,
therefore we assume $\alpha$ = 0 and $K(\alpha)$ = 1. The estimated exchange
interactions for each Gd3+ concentration and applied magnetic field are
summarized in Table 1. Although we should underestimate $J_{fp}(0)$ and
$<J_{fp}^{2}(\textbf{q})>^{1/2}$ due to remaining CEF effects and a local
reduction of the DOS, the trend observed in our analysis qualitatively
trustworthy.

##  IV IV. Discussion

As clearly shown in Table 1, in X-band measurements we do have a clear
q-dependence in the exchange interaction for (Bi1-xGdx)2Se3, which appears to
reduce upon increasing the applied magnetic field (Q-band). A notion of this
effect can be obtained looking to the real space. Our results indicate that,
at lower fields, there is a stronger exchange interaction surrounding the Gd3+
ions. In other words, it appears that occurs a localization of the carriers
surrounding the Gd3+ ions. The increase of the magnetic field appears to
induce a magnetic breakdown to the effect, and the system starts to behave
more like a regular metal with a simple Fermi surface with no q-dependence.
During this whole process, we should expect that the average of the
interaction between conduction electrons and 4$f$ local moments is almost
constant, which explains why there is an evolution of $\Delta g$ while $b$
remains unchanged. The last important piece of information to pinpoint our
interpretation comes from the Gd3+ spectra evolution as a function of the
applied magnetic field. The transition from an insulating-like Gd3+ spectrum
at X-band to a collapse of the CEF splitting at Q-band data indicates a
destructive-like interference of the carriers surrounding the Gd3+ ions, which
is the so-called local WAL effect.

A direct consequence of the local WAL effect, we may also understand the
collapse of the CEF splitting of the Gd3+ ions from a viewpoint of CEF
contributions from the lattice and the carriers. The CEF parameters may be
affected by the surrounding charges at the Gd3+ site and eventually they can
be strongly reduced. In other words, it appears that, due to the $S$-ion
character of Gd3+ ions, the influence of the conduction electrons could become
relevant to the CEF effects compared to the lattice charges, as already shown
in half-Heusler systems [86, 87]. As such, a change in the local charge
distribution could cause a collapse of the Gd3+ fine structure due to a strong
reduction of the crystal field parameters associated with two contributions of
different signs of the CEF parameters (lattice and carriers contributions).
Since the second order crystal field parameter $b_{2}^{0}$ is much larger than
the fourth order ones in (Bi1-xGdx)2Se3 [47], we can estimate $|b_{2}^{0}|$ =
32/3*12 $\leq$ 0.9 Oe for Q-band measurements [88].

One important concern about this proposed scenario is the comparison with
previous transport results. Albeit previous reports show that even bulk Bi2Se3
samples have a WAL effect at low fields, and thinner samples have a more
significant critical field ($H$ $\leq$ 1 T) [50], we must pay attention that
all the results report WAL effects at low temperatures. In order to understand
the origin of this discrepancy, it is important to understand that transport
and ESR are different techniques. While the response in transport is a
macroscopic, global property, in ESR we obtain a microscopic viewpoint. In ESR
measurements we have two distinct relaxation channels: the relaxation through
the spin-phonon process and the spin-spin relaxation, which is faster and
involves the carriers of the system [61, 62, 63]. From our Gd3+ spin dynamics,
Fig. 3, we clearly obtain signatures that the coupling between the 4$f$ local
moments and the carriers are relevant, and they should dominate the relaxation
of the system. Concomitantly, Gd3+ has $L$ = 0, which means that the spin-
phonon coupling is rather weak. That means that any potential scattering
involving phonons, which would mask the WAL signatures, are going to be
heavily suppressed. Such locality also help us to have an insight of the
difference of critical fields. For thicker transport samples, the cusp has a
critical field of only $\approx$ 1 kOe, while our X-band measurements occur at
fields of 3.5 kOe. Again, the phonons contribution are going to be heavily
suppressed, which means that it is understandable that the critical field for
the local WAL effect might be comparable to the fields of thinner samples.
Nonetheless, it is worth to note that Q-band measurements have fields applied
at the order of $\approx$ 12 kOe. An alternate scenario could also rely on
strain effects, which are originated from surface effects, playing a bigger
role in Q-band measurements. In this scenario, the central line could have a
smaller line width and we would be able to describe the data without an
additional collapsed spectrum. Two different results indicate that this is an
unlikely explanation. First of all, the data is better described with two
resonances of different line widths at the same $g$-value. Even with we do not
take this fact into account, the intensity of the central line should be, at
least, double than the expected value to fairly describe the data - which is
inconsistent with the crystal field Hamiltonian. Another point is that the
skin depth is $\sim$ 11 $\mu$m and 8 $\mu$m for X- and Q-band measurements,
respectively. With those skin depths, as expected, we do not see any
signatures of surface effects in our data. Therefore, as far as our
experimental data shows, the bulk dominates the ESR signatures and strain
effects should not play a role.

It is worth also pointing the evolution of the sign of the $g$-shift. As
already mentioned, our results indicate that $p$-type carriers are the main
contributors near the Fermi level [73, 74]. However, earlier reports show
positive $g$-shifts (typical for $s$-type carriers) for Gd3+ substituted
Bi2(Te,Se)3 [46, 47]. Perhaps, the higher temperature reached in their
synthesis increases the number of defects, affecting the $g$-value - which is
expected due to the Se and Te boiling points. Future systematic studies of
Gd3+ substituted Bi2Se3 synthesized at different temperatures could help to
enlighten this open question.

The defects, mainly Se vacancies, also appears to have an important role in
the interpretation of our experiments. First of all, we still observed a weak,
but visible, Gd3+ fine structure at Q-band measurements. Due to the different
CEF states, Gd3+ local moments closer to impurity centers may need higher
fields to completely suppress the local WAL effect. Albeit it should be taken
with care, it is possible to estimate the ratio between Gd3+ ions in an
insulating and conducting environment performing a double integration from our
Gd3+ ESR spectrum for Q-band measurements. We roughly estimate that 70 % of
the Gd3+ ions have a conducting-like environment. As mentioned before,
previous results indicate that vacancies tend to be inhomogeneous in Bi2Se3
matrix [71]. Although the distribution of Gd3+ ions appears to be homogeneous,
it appears that those domains may have a tendency of accumulating near the
magnetic impurities, which would be consistent with a high number of local
moments center still showing an insulating character. The second hint of the
role of vacancies comes from the Gd3+ spin dynamics, especially the changes
observed in the ESR data at low temperatures when comparing X-band and Q-band
data. In this temperature range, presumably, the interaction between impurity
centers and Gd3+ local moments is more relevant, We obtained an evolution of
the negative increase of the $g$-shift at low temperatures as a function of
field: the Gd3+ $g$-shift decrease due to spin-spin interactions is much more
pronounced at X-band measurements when compared to Q-band ones [Fig. 3 c)].
The change in the charge distribution would change the interaction between
defects and 4$f$ local moments and, naturally, we obtain a change in the
low-$T$ Gd3+ spin dynamics.

Regarding the topology of the system, although the strong spin-orbit coupling
of the bulk is essential to the WAL mechanism, it is not clear if the topology
of the system plays any role. Another potential signature of surface
excitations could appear in a diffusive-like contribution to the Gd3+ ESR line
shape [89, 90, 91]. However, the presence of carriers will naturally suppress
any contribution from the surface excitations to the relaxation. This helps to
understand, associated with the weak CEF effects of Gd3+ ions, the lack of any
diffusive-like effects in the Gd3+ line shape [89, 90, 91].

Going back to the magnetism of 4$f$ electrons in chalcogenides, the obtained
exchange interactions, which is a Ruderman-Kittel-Kasuya-Yosida (RKKY)
interaction, in Table 1 are clearly smaller than those obtained for
Bi2-xMnxTe3 [23]. Additionally, it has been observed a change of the Korringa
rate as a function of Mn2+ concentration, which may be linked to changes in
the $p-d$ hybridization. In the Gd3+ case, we observed the same Korringa
relaxation for all the explored concentrations, indicating that $f-p$
hybridization may only play a small role, highlighting the clear difference
between the effects of $d$ and $f$ substitutions.

##  V V. Conclusion

In summary, we performed electron spin resonance and complementary macroscopic
measurements in the Gd3+-substituted Bi2Se3 ($x$ = 0, 0.001, 0.002 and 0.006).
The Gd3+ ESR spectra at 4 K for different concentrations show seven Dysonian
lines for X-band measurements, which evolves to an apparent contribution from
distinct Gd3+ sites spectra for Q-band. We conjecture that such evolution
might be due to a breakdown of the local WAL effect and a change in the
crystal field parameters in the vicinity of the vacancies. This interpretation
is consistent with the Gd3+ spin dynamics response. Additionally, we show that
the 4$f$ substitution does not increase the DOS at the Fermi level, neither
introduces a relevant $f-p$ hybridization. This is in contrast to more
traditional substitutions with transition metals $d$ states magnetic ions and
indicates different magnetic mechanisms for $d$ and $f$ states in this system.
Our work points out that 4$f$ substitution in chalcogenides is an interesting
path to explore even further the role of magnetic impurities in this model
system.

###### Acknowledgements.

We would like to thank Dr. S. K. Misra for providing Ref. [88]. This work was
supported by FAPESP (SP-Brazil) Grants No 2022/09240-3, 2020/12283-0,
2019/26247-9, 2018/11364-7, 2017/10581-1, 2012/04870-7, CNPq Grants No
311783/2021-0, 309483/2018-2, 442230/2014-1 and 304649/2013-9, CAPES and
FINEP-Brazil. Work at Los Alamos was supported by the Los Alamos Laboratory
Directed Research and Development program through project 20210064DR.

## References

  * Bernevig  _et al._ [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006).
  * König  _et al._ [2007] M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007).
  * Culcer  _et al._ [2020] D. Culcer, A. C. Keser, Y. Li, and G. Tkachov, 2D Mater. 7, 022007 (2020).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Tokura  _et al._ [2017] Y. Tokura, M. Kawasaki, and N. Nagaosa, Nat. Phys. 13, 1056 (2017).
  * Zhang  _et al._ [2012] D. Zhang, A. Richardella, D. W. Rench, S.-Y. Xu, A. Kandala, T. C. Flanagan, H. Beidenkopf, A. L. Yeats, B. B. Buckley, P. V. Klimov, _et al._ , Phys. Rev. B 86, 205127 (2012).
  * Chang  _et al._ [2013] C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, _et al._ , Science 340, 167 (2013).
  * He  _et al._ [2014] K. He, Y. Wang, and Q.-K. Xue, Natl. Sci. Rev. 1, 38 (2014).
  * Chang  _et al._ [2015] C.-Z. Chang, W. Zhao, D. Y. Kim, H. Zhang, B. A. Assaf, D. Heiman, S.-C. Zhang, C. Liu, M. H. Chan, and J. S. Moodera, Nat. Mater. 14, 473 (2015).
  * Grauer  _et al._ [2015] S. Grauer, S. Schreyeck, M. Winnerlein, K. Brunner, C. Gould, and L. Molenkamp, Phys. Rev. B 92, 201304 (2015).
  * Xiao  _et al._ [2018] D. Xiao, J. Jiang, J.-H. Shin, W. Wang, F. Wang, Y.-F. Zhao, C. Liu, W. Wu, M. H. Chan, N. Samarth, _et al._ , Phys. Rev. Lett. 120, 056801 (2018).
  * Yue  _et al._ [2019] C. Yue, Y. Xu, Z. Song, H. Weng, Y.-M. Lu, C. Fang, and X. Dai, Nat. Phys. 15, 577 (2019).
  * Zhang  _et al._ [2019] D. Zhang, M. Shi, T. Zhu, D. Xing, H. Zhang, and J. Wang, Phys. Rev. Lett. 122, 206401 (2019).
  * Deng  _et al._ [2020] Y. Deng, Y. Yu, M. Z. Shi, Z. Guo, Z. Xu, J. Wang, X. H. Chen, and Y. Zhang, Science 367, 895 (2020).
  * Liu  _et al._ [2020] C. Liu, Y. Wang, H. Li, Y. Wu, Y. Li, J. Li, K. He, Y. Xu, J. Zhang, and Y. Wang, Nat. mater. 19, 522 (2020).
  * Nenno  _et al._ [2020] D. M. Nenno, C. A. Garcia, J. Gooth, C. Felser, and P. Narang, Nat. Rev. Phys. 2, 682 (2020).
  * Fijalkowski  _et al._ [2021a] K. Fijalkowski, N. Liu, M. Hartl, M. Winnerlein, P. Mandal, A. Coschizza, A. Fothergill, S. Grauer, S. Schreyeck, K. Brunner, _et al._ , Phys. Rev. B 103, 235111 (2021a).
  * Zhang  _et al._ [2009] H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Nature physics 5, 438 (2009).
  * Analytis  _et al._ [2010] J. G. Analytis, J.-H. Chu, Y. Chen, F. Corredor, R. D. McDonald, Z. Shen, and I. R. Fisher, Physical Review B 81, 205407 (2010).
  * Pan  _et al._ [2011] Z.-H. Pan, E. Vescovo, A. Fedorov, D. Gardner, Y. Lee, S. Chu, G. Gu, and T. Valla, Physical review letters 106, 257004 (2011).
  * Fijalkowski  _et al._ [2021b] K. M. Fijalkowski, N. Liu, P. Mandal, S. Schreyeck, K. Brunner, C. Gould, and L. W. Molenkamp, Nature communications 12, 1 (2021b).
  * Mahani  _et al._ [2014] M. R. Mahani, A. Pertsova, M. F. Islam, and C. M. Canali, Phys. Rev. B 90, 195441 (2014).
  * Zimmermann  _et al._ [2016] S. Zimmermann, F. Steckel, C. Hess, H. Ji, Y. S. Hor, R. J. Cava, B. Büchner, and V. Kataev, Phys. Rev. B 94, 125205 (2016).
  * Peixoto  _et al._ [2016] T. R. Peixoto, H. Bentmann, S. Schreyeck, M. Winnerlein, C. Seibel, H. Maaß, M. Al-Baidhani, K. Treiber, S. Schatz, S. Grauer, _et al._ , Phys. Rev. B 94, 195140 (2016).
  * Zhang  _et al._ [2018] W. Zhang, D. West, S. H. Lee, Y. Qiu, C.-Z. Chang, J. S. Moodera, Y. San Hor, S. Zhang, and W. Wu, Phys. Rev. B 98, 115165 (2018).
  * Ye  _et al._ [2019] M. Ye, T. Xu, G. Li, S. Qiao, Y. Takeda, Y. Saitoh, S.-Y. Zhu, M. Nurmamat, K. Sumida, Y. Ishida, _et al._ , Phys. Rev. B 99, 144413 (2019).
  * Bouaziz  _et al._ [2019] J. Bouaziz, M. dos Santos Dias, F. S. M. Guimarães, and S. Lounis, Phys. Rev. Mater. 3, 054201 (2019).
  * Tcakaev  _et al._ [2020a] A. Tcakaev, V. Zabolotnyy, R. Green, T. Peixoto, F. Stier, M. Dettbarn, S. Schreyeck, M. Winnerlein, R. C. Vidal, S. Schatz, _et al._ , Phys. Rev. B 101, 045127 (2020a).
  * Peixoto  _et al._ [2020] T. R. Peixoto, H. Bentmann, P. Rüßmann, A.-V. Tcakaev, M. Winnerlein, S. Schreyeck, S. Schatz, R. C. Vidal, F. Stier, V. Zabolotnyy, _et al._ , npj Quantum Mater. 5, 1 (2020).
  * Islam  _et al._ [2018] M. Islam, C. M. Canali, A. Pertsova, A. Balatsky, S. Mahatha, C. Carbone, A. Barla, K. Kokh, O. Tereshchenko, E. Jiménez, _et al._ , Phys. Rev. B 97, 155429 (2018).
  * Mogi  _et al._ [2017a] M. Mogi, M. Kawamura, R. Yoshimi, A. Tsukazaki, Y. Kozuka, N. Shirakawa, K. Takahashi, M. Kawasaki, and Y. Tokura, Nat. Mater. 16, 516 (2017a).
  * Mogi  _et al._ [2017b] M. Mogi, M. Kawamura, A. Tsukazaki, R. Yoshimi, K. S. Takahashi, M. Kawasaki, and Y. Tokura, Sci. Adv. 3, eaao1669 (2017b).
  * Mathimalar  _et al._ [2020] S. Mathimalar, S. Sasmal, A. Bhardwaj, S. Abhaya, R. Pothala, S. Chaudhary, B. Satpati, and K. V. Raman, npj Quantum Mater. 5, 1 (2020).
  * Pereira  _et al._ [2020a] V. Pereira, S. Altendorf, C. Liu, S. Liao, A. Komarek, M. Guo, H.-J. Lin, C. Chen, M. Hong, J. Kwo, _et al._ , Phys. Rev. Mater. 4, 064202 (2020a).
  * Pereira  _et al._ [2020b] V. Pereira, C.-N. Wu, C.-A. Knight, A. Choa, L. Tjeng, and S. Altendorf, APL Mater. 8, 071114 (2020b).
  * Chen  _et al._ [2015] T. Chen, W. Liu, F. Zheng, M. Gao, X. Pan, G. Van Der Laan, X. Wang, Q. Zhang, F. Song, B. Wang, _et al._ , Adv. Mater. 27, 4823 (2015).
  * Tcakaev  _et al._ [2020b] A. Tcakaev, V. B. Zabolotnyy, C. I. Fornari, P. Rüßmann, T. R. Peixoto, F. Stier, M. Dettbarn, P. Kagerer, E. Weschke, E. Schierle, _et al._ , Phys. Rev. B 102, 184401 (2020b).
  * Li  _et al._ [2013] S. Li, S. Harrison, Y. Huo, A. Pushp, H. Yuan, B. Zhou, A. Kellock, S. Parkin, Y.-L. Chen, T. Hesjedal, _et al._ , Applied Physics Letters 102, 242412 (2013).
  * Filnov  _et al._ [2020] S. Filnov, I. I. Klimovskikh, D. Estyunin, A. Fedorov, V. Y. Voroshnin, A. Koroleva, A. G. Rybkin, E. Shevchenko, Z. S. Aliev, M. Babanly, _et al._ , Phys. Rev. B 102, 085149 (2020).
  * Lee  _et al._ [2019] E. Lee, S. Seong, M. Y. Yang, J. Kim, M.-H. Jung, B.-G. Park, Y. Kim, S. W. Han, and J.-S. Kang, Appl. Phys. Lett. 115, 072404 (2019).
  * Lee  _et al._ [2015] H. Lee, J. Kim, K. Lee, A. Jelen, S. Vrtnik, Z. Jagličić, J. Dolinšek, and M. Jung, Appl. Phys. Lett. 107, 182409 (2015).
  * Kim and Jung [2018] S.-W. Kim and M.-H. Jung, AIP Adv. 8, 101319 (2018).
  * Filnov  _et al._ [2019] S. Filnov, Y. A. Surnin, A. Koroleva, I. Klimovskikh, D. Estyunin, A. Y. Varykhalov, K. Bokai, K. Kokh, O. Tereshchenko, V. Golyashov, _et al._ , J. Exp. Theor. Phys. 129, 404 (2019).
  * Kholdi  _et al._ [1994] M. E. Kholdi, M. Averous, S. Charar, C. Fau, G. Brun, H. Ghoumari-Bouanani, and J. Deportes, Phys. Rev. B 49, 1711 (1994).
  * Song  _et al._ [2012] Y. Song, F. Yang, M.-Y. Yao, F. Zhu, L. Miao, J.-P. Xu, M.-X. Wang, H. Li, X. Yao, F. Ji, _et al._ , Appl. Phys. Lett. 100, 242403 (2012).
  * Isber  _et al._ [1995] S. Isber, S. Charar, V. Mathet, C. Fau, and M. Averous, Phys. Rev. B 51, 15578 (1995).
  * Gratens  _et al._ [1997] X. Gratens, S. Isber, S. Charar, C. Fau, M. Averous, S. K. Misra, Z. Golacki, M. Ferhat, and J. C. Tedenac, Physical Review B 55, 8075 (1997).
  * Garitezi  _et al._ [2015] T. Garitezi, G. Lesseux, C. Jesus, T. Grant, Z. Fisk, R. Urbano, C. Rettori, and P. Pagliuso, in  _Journal of Physics: Conference Series_ , Vol. 592 (IOP Publishing, 2015) p. 012125.
  * Hikami  _et al._ [1980] S. Hikami, A. I. Larkin, and Y. Nagaoka, Prog. Theor. Phys. 63, 707 (1980).
  * Kim  _et al._ [2011] Y. S. Kim, M. Brahlek, N. Bansal, E. Edrey, G. A. Kapilevich, K. Iida, M. Tanimura, Y. Horibe, S.-W. Cheong, and S. Oh, Phys. Rev. B 84, 073109 (2011).
  * Lu and Shen [2011] H.-Z. Lu and S.-Q. Shen, Phys. Rev. B 84, 125138 (2011).
  * Lu  _et al._ [2011] H.-Z. Lu, J. Shi, and S.-Q. Shen, Phys. Rev. Lett. 107, 076801 (2011).
  * Huang  _et al._ [2012] F.-T. Huang, M.-W. Chu, H. Kung, W. Lee, R. Sankar, S.-C. Liou, K. Wu, Y. Kuo, and F. Chou, Phys. Rev. B 86, 081104 (2012).
  * Devidas  _et al._ [2015] T. Devidas, E. Amaladass, S. Sharma, R. Rajaraman, D. Sornadurai, N. Subramanian, A. Mani, C. Sundar, and A. Bharathi, EPL (Europhys. Lett.) 108, 67008 (2015).
  * Adroguer  _et al._ [2015] P. Adroguer, W. E. Liu, D. Culcer, and E. Hankiewicz, Phys. Rev. B 92, 241402 (2015).
  * Teng  _et al._ [2019] J. Teng, N. Liu, and Y. Li, J. Semicond. 40, 081507 (2019).
  * Janíček  _et al._ [2008] P. Janíček, Č. Drašar, P. Lošták, J. Vejpravová, and V. Sechovskỳ, Phys. B: Condens. Matter 403, 3553 (2008).
  * Chen  _et al._ [2011] J. Chen, X. He, K. Wu, Z. Ji, L. Lu, J. Shi, J. Smet, and Y. Li, Phys. Rev. B 83, 241304 (2011).
  * Wang  _et al._ [2011] J. Wang, A. M. DaSilva, C.-Z. Chang, K. He, J. Jain, N. Samarth, X.-C. Ma, Q.-K. Xue, and M. H. Chan, Phys. Rev. B 83, 245438 (2011).
  * Pal  _et al._ [2012] H. Pal, V. Yudson, and D. Maslov, Phys. Rev. B 85, 085439 (2012).
  * Abragam and Bleaney [2012] A. Abragam and B. Bleaney,  _Electron paramagnetic resonance of transition ions_ (OUP Oxford, 2012).
  * Barnes [1981] S. Barnes, Adv. Phys. 30, 801 (1981).
  * Poole and Farach [1971] C. P. Poole and H. A. Farach,  _Relaxation in magnetic resonance_ , Vol. 19 (Elsevier, 1971).
  * Feher and Kip [1955] G. Feher and A. Kip, Phys. Rev. 98, 337 (1955).
  * Dyson [1955] F. J. Dyson, Phys. Rev. 98, 349 (1955).
  * Anderson and Weiss [1953] P.-W. Anderson and P. Weiss, Reviews of Modern Physics 25, 269 (1953).
  * Urban  _et al._ [1975] P. Urban, D. Davidov, B. Elschner, T. Plefka, and G. Sperlich, Phys. Rev. B 12, 72 (1975).
  * Park  _et al._ [2013] K. Park, Y. Nomura, R. Arita, A. Llobet, and D. Louca, Phys. Rev. B 88, 224108 (2013).
  * Xu  _et al._ [2014] C. Xu, A. Hewitt, J. Wang, T. Guan, J. Boltersdorf, P. A. Maggard, D. B. Dougherty, and K. Gundogdu, J. Appl. Phys. 116, 043519 (2014).
  * Tayal  _et al._ [2017] A. Tayal, D. Kumar, and A. Lakhani, J. Phys.: Condens. Matter 29, 445704 (2017).
  * Taylor  _et al._ [2012] R. E. Taylor, B. Leung, M. P. Lake, and L.-S. Bouchard, J. Phys. Chem. C 116, 17300 (2012).
  * Davidov  _et al._ [1973] D. Davidov, C. Rettori, A. Dixon, K. Baberschke, E. Chock, and R. Orbach, Phys. Rev. B 8, 3563 (1973).
  * Cao  _et al._ [2013] Y. Cao, J. Waugh, X. Zhang, J.-W. Luo, Q. Wang, T. Reber, S. Mo, Z. Xu, A. Yang, J. Schneeloch, _et al._ , Nat. Phys. 9, 499 (2013).
  * Vidal  _et al._ [2013] F. Vidal, M. Eddrief, B. R. Salles, I. Vobornik, E. Velez-Fort, G. Panaccione, and M. Marangolo, Physical Review B 88, 241410 (2013).
  * Kim  _et al._ [2015] S. Kim, S. Vrtnik, J. Dolinšek, and M. Jung, Appl. Phys. Lett. 106, 252401 (2015).
  * Shimizu  _et al._ [2006] Y. Shimizu, K. Miyagawa, K. Kanoda, M. Maesato, and G. Saito, Phys. Rev. B 73, 140407 (2006).
  * Riedl  _et al._ [2019] K. Riedl, R. Valentí, and S. M. Winter, Nat. Commun. 10, 1 (2019).
  * Kawamura and Uematsu [2019] H. Kawamura and K. Uematsu, J. Phys. Condens. Matter 31, 504003 (2019).
  * Pustogow  _et al._ [2020] A. Pustogow, T. Le, H.-H. Wang, Y. Luo, E. Gati, H. Schubert, M. Lang, and S. Brown, Phys. Rev. B 101, 140401 (2020).
  * Miksch  _et al._ [2021] B. Miksch, A. Pustogow, M. J. Rahim, A. A. Bardin, K. Kanoda, J. A. Schlueter, R. Hübner, M. Scheffler, and M. Dressel, Science 372, 276 (2021).
  * Rettori  _et al._ [1974] C. Rettori, H. Kim, E. Chock, and D. Davidov, Phys. Rev. B 10, 1826 (1974).
  * Moriya [1963] T. Moriya, Journal of the Physical Society of Japan 18, 516 (1963).
  * Narath [1967] A. Narath, Phys. Rev. 163, 232 (1967).
  * Narath and Weaver [1968] A. Narath and H. Weaver, Phys. Rev. 175, 373 (1968).
  * Shaw Jr and Warren Jr [1971] R. W. Shaw Jr and W. W. Warren Jr, Phys. Rev. B 3, 1562 (1971).
  * Pagliuso  _et al._ [1999] P. Pagliuso, C. Rettori, M. Torelli, G. Martins, Z. Fisk, J. Sarrao, M. Hundley, and S. Oseroff, Phys. Rev. B 60, 4176 (1999).
  * Souza  _et al._ [2019] J. Souza, C. Jesus, G. Lesseux, P. Rosa, R. Urbano, and P. Pagliuso, J. Phys. Condens. Matter 31, 465701 (2019).
  * Misra [1986] S. K. Misra, Magn. Reson. Rev. 10, 285 (1986).
  * Lesseux  _et al._ [2016] G. Lesseux, T. Garitezi, P. Rosa, C. Jesus, S. Oseroff, J. Sarrao, Z. Fisk, R. Urbano, P. Pagliuso, and C. Rettori, J. Phys. Condens. Matter 28, 125601 (2016).
  * Souza  _et al._ [2018] J. Souza, G. Lesseux, R. Urbano, C. Rettori, and P. Pagliuso, AIP adv. 8, 055713 (2018).
  * Souza  _et al._ [2021] J. Souza, M. König, M. Ale Crivillero, M. Malcolms, R. Urbano, Z. Fisk, P. Rosa, P. Pagliuso, S. Wirth, and J. Sichelschmidt, Phys. Rev. Research 3, 033016 (2021).