# Interplay between breathing-mode distortions and magnetic order in rare-
earth nickelates from ab initio magnetic models

Danis I. Badrtdinov  Theoretical Physics and Applied Mathematics Department,
Ural Federal University, 620002 Yekaterinburg, Russia  Center for
Computational Quantum Physics, Flatiron Institute, 162 Fifth avenue, New York,
NY 10010, USA  Alexander Hampel<EMAIL_ADDRESS>Center for
Computational Quantum Physics, Flatiron Institute, 162 Fifth avenue, New York,
NY 10010, USA  Cyrus E. Dreyer  Center for Computational Quantum Physics,
Flatiron Institute, 162 Fifth avenue, New York, NY 10010, USA  Department of
Physics and Astronomy, Stony Brook University, Stony Brook, New York,
11794-3800, USA

###### Abstract

We use density-functional theory calculations to explore the magnetic
properties of perovskite rare-earth nickelates, $\mathcal{R}$NiO3, by
constructing microscopic magnetic models containing all relevant exchange
interactions via Wannierization and Green’s function techniques. These models
elucidate the mechanism behind the formation of antiferromagnetic order with
the experimentally observed propagation vector, and explain the reason
previous DFT plus Hubbard $U$ calculations favored ferromagnetic order. We
perform calculations of magnetic moments and exchange-coupling parameters for
different amplitudes of the $R_{1}^{+}$ breathing mode distortion, which
results in expanded and compressed NiO6 octahedra. We find that the magnetic
moment vanishes for the “short bond” nickels, i.e., the ones in the compressed
octahedra. The inclusion of spin-orbit coupling demonstrates that the magnetic
anisotropy is very small, while the magnetic moment of the short bond nickel
atoms tend to zero even for the noncollinear case. Our results provide a clear
picture of the trends of the magnetic order across the nickelate series and
give insights into the coupling between magnetic order and structural
distortions.

##  I Introduction

Rare-earth (RE) nickelate perovskites, with the chemical formula
$\mathcal{R}$NiO3, have received significant attention due to their rich phase
diagram. As a function of the ionic size of RE element $\mathcal{R}$ Hampel
and Ederer (2017); Varignon  _et al._ (2017); Catalano  _et al._ (2018); Lu
_et al._ (2017), a combined structural and metal-insulator transition (MIT)
occurs with temperature, as well as the formation of long-range magnetic
order. Both can be tuned by external pressure Obradors  _et al._ (1993),
strain Catalano  _et al._ (2014) and reduced dimensionality Boris  _et al._
(2011). This tunability paves the way for engineering nanoelectronic and
spintronic devices Y  _et al._ (2015); Shi  _et al._ (2014); Cheong and
Mostovoy (2007).

The structural transition, which occurs in all RE nickelates except LaNiO3, is
characterized by a symmetry lowering from the orthorhombic $Pbnm$ to the
monoclinic $P2_{1}/n$ spacegroup upon cooling; it can be described mainly by a
“breathing-mode” distortion, which creates expanded and compressed NiO6
octahedra arranged in a three dimensional checkerboard pattern. We refer to
the two different Ni sites that result as residing in either the long bond
(LB) or short bond (SB) octahedra Alonso  _et al._ (1999) (Fig. 1).

This distortion strongly modifies the electronic configuration and
corresponding magnetic moments of nickel atoms. One of the proposed
descriptions of the resulting electronic structure Mazin  _et al._ (2007)
suggests that the initial electronic configuration of the Ni atoms, i.e.,
$3d^{7}$ with spin $S=\frac{1}{2}$, transforms to two charge-disproportionated
Ni sublattices $3d^{6}\,(S=0)$ and $3d^{8}\,(S=1)$ under the distortion.
However, this atomic picture does not consider the strong hybridization
between Ni-$3d$ and O-$2p$ orbitals. Hence, later studies introduced a more
complete description of the electronic configuration including oxygen (ligand)
“holes,” denoted as $\underline{L}$ Johnston  _et al._ (2014); Green  _et al._
(2016), and describing the connected structural-MIT in nickelates as a site-
selective Mott transition Park  _et al._ (2012); Lau and Millis (2013);
Johnston  _et al._ (2014); Subedi  _et al._ (2015); Haule and Pascut (2017);
Peil  _et al._ (2019). The resulting SB and LB octahedra are described by
$d^{8}\underline{L}^{0}$ and $d^{8}\underline{L}^{2}$ states with total spins
tending toward $S=0$ and $S=1$, respectively. Conventional density functional
theory (DFT) based methods also support this interpretation, showing different
magnetic moments of NiSB and NiLB atoms, while demonstrating nearly equal
integrated charges around both Ni atoms inside the projector augmented wave
(PAW) spheres Park  _et al._ (2012); Hampel and Ederer (2017). However, such
an analysis will depend on the size of the augmentation region Kresse and
Joubert (1999), and cannot fully capture the hybridization with the oxygens; a
more complete picture may be obtained, e.g., via constructing Wannier
functions for the (possibly hybridized) Ni $3d$ states Varignon  _et al._
(2017).

At low temperatures, RE nickelates (again with the exception of LaNiO3)
exhibit an antiferromagnetic (AFM) insulating ground state. The magnetic
structure has a propagation vector of
$\mathbf{q}=(\frac{1}{4},\frac{1}{4},\frac{1}{4})$ (in pseudocubic notation),
which was probed by neutron scattering for bulk systems García-Muñoz  _et al._
(1994); Alonso  _et al._ (1999); Gawryluk  _et al._ (2019); Muñoz  _et al._
(2009); Fernández-Díaz  _et al._ (2001) and resonant soft X-ray scattering
(RIXS) for thin film heterostructures Scagnoli  _et al._ (2006); Lu  _et al._
(2018). The exact arrangement of magnetic moments in the systems is still
under debate and strongly depends on the interpretation of experimental data
Haule and Pascut (2017); Catalano  _et al._ (2018). Particularly, neutron data
suggests a collinear magnetic ordering, where the magnetic moment of SB nickel
atoms vanish Muñoz  _et al._ (2009); Fernández-Díaz  _et al._ (2001); Haule
and Pascut (2017). One candidate for the ground-state magnetic configuration
is a T-type AFM (T-AFM) structure Catalano  _et al._ (2018); Hampel and Ederer
(2017), depicted in Fig. 2(a). In this structure, the LB moments are arranged
in an “up-up-down-down” ($\uparrow\uparrow\downarrow\downarrow$) pattern in
all three Cartesian directions, and the SB moments tend to zero with
increasing $R_{1}^{+}$ amplitude.

On the other hand, the corresponding RIXS spectra have been interpreted as a
non-collinear (NCL) arrangement of spins with clearly non-zero SB magnetic
moments Scagnoli  _et al._ (2006); Lu  _et al._ (2018), also involving the
$\uparrow\uparrow\downarrow\downarrow$ pattern in all three Cartesian
directions [see Fig. 2(b)]. Here, all moments lie in the $ac$ plane, with
their direction reversed upon moving along [1, 0, 1]. This magnetic
configuration can be considered as an orthogonal spin spiral. More
importantly, the nearest neighbor SB and LB magnetic moments are perpendicular
to each other, which can weaken interactions between them.

Our study is aimed at understanding the formation of these proposed complex
magnetic orderings on a microscopic level in the RE nickelates. To this end we
performed a systematic DFT plus Hubbard $U$ (DFT+$U$) analysis of the magnetic
properties of a series of RE nickelates: PrNiO3, SmNiO3 and LuNiO3. We will
consider the collinear T-AFM and NCL order [Fig. 2(a) and (b)], as well as the
ferromagnetic (FM) case, which was shown to have the lowest energy in previous
DFT+$U$ studies Hampel and Ederer (2017); Varignon  _et al._ (2017). Spin
models are constructed by obtaining Wannier functions for the hybridized
Ni($3d$)-O($2p$) orbitals, and then using the local force theorem approach
Liechtenstein  _et al._ (1987); Mazurenko and Anisimov (2005) to extract
exchange couplings. We analyze the evolution of the magnetic interactions
under the $R_{1}^{+}$ structural breathing mode distortion. Our calculations
reveal strong ferromagnetic nearest-neighbor and weaker antiferromagnetic
next-nearest-neighbor couplings, which are induced due to magnetically active
Ni $d_{z^{2}}$ and $d_{x^{2}-y^{2}}$ states. Competition between these
couplings plays the dominant role in determining the formation of magnetic
order in $\mathcal{R}$NiO3 systems and explains why previous theoretical
calculations favored ferromagnetic ordering. The inclusion of spin-orbit
coupling demonstrates that the magnetic anisotropy is very small, while the
magnetic moments of SB nickel atoms tend to zero for finite $R_{1}^{+}$
amplitude even in the NCL case.

The remainder of this paper is organized as follows. Sec. II covers the
methods and theories used in our calculations. In Sec. III, we present our
main results for magnetic moments and interactions of selected rare-earth
nickelate systems: PrNiO3, SmNiO3 and LuNiO3. We discuss the implications of
our results on the proposed magnetic structures in Sec. IV. Finally, Sec. V
concludes the paper.

Figure 1:  Evolution of crystal structure of $\mathcal{R}$NiO3 under the
breathing mode distortion $R_{1}^{+}$, which leads to formation of short bond
(SB) and long bond (LB) NiO6 octahedra. The rare-earth element $\mathcal{R}$
is denoted by yellow spheres. Crystal structures are visualized using the
VESTA software Momma and Izumi (2011).

##  II Theoretical framework

###  II.1 Formalism

We construct the following spin model, which describes energetics of
individual spins ${\bf S}_{i}$ on nickel atoms interacting via isotropic
exchange couplings $J_{ij}$:

$\hat{\mathcal{H}}^{\text{spin}}=\sum_{i>j}J_{ij}\hat{{\bf S}}_{i}\hat{{\bf
S}}_{j}.$ (1)

The exchange integrals can be calculated using the local force theorem
approach Liechtenstein  _et al._ (1987); Mazurenko and Anisimov (2005)

$\begin{split}J_{ij}&=\frac{1}{2\pi
S^{2}}\sum\limits_{m,m^{\prime},n,n^{\prime}}\times\\\ &{\rm
Im}\int\limits_{-\infty}^{E_{F}}d\epsilon\left[\Delta^{mm^{\prime}}_{i}G^{m^{\prime}n}_{ij\downarrow}(\epsilon)\Delta^{nn^{\prime}}_{j}G^{n^{\prime}m}_{ji\uparrow}(\epsilon)\right],\end{split}$
(2)

where $m,m^{\prime},n,n^{\prime}$ are orbital indices running over Ni 3$d$
$(l=2)$ states, $E_{F}$ is the Fermi energy, $S=\frac{1}{2}$ is the spin
quantum number 111For simplicity, we use spin quantum number $S=\frac{1}{2}$
in Eq. (2) during the whole range of breathing distortions. It was revealed
that very large breathing distortion amplitudes are needed to achieve spins
$0$ and $1$ for SB and LB magnetic moments Johnston  _et al._ (2014); Green
_et al._ (2016), which is much higher than our considered cases,
$\Delta^{mm^{\prime}}_{i}=H^{mm^{\prime}}_{ii\uparrow}-H^{mm^{\prime}}_{ii\downarrow}$
is the on-site potential, and $G$ is the single-particle Green’s functions

$\displaystyle
G^{nn^{\prime}}_{ij\sigma}(\epsilon)=\frac{1}{N_{\textbf{k}}}\sum\limits_{\mathbf{k},l}^{M}\frac{c^{nl}_{i\sigma}(\mathbf{k})c^{n^{\prime}l^{*}}_{j\sigma}(\bf{k})}{\epsilon-E^{l}_{\sigma}(\mathbf{k})}.$
(3)

Here $N_{\textbf{k}}$ is the number of $k$ points in the first Brillouin zone
(BZ), $c^{n^{\prime}l*}_{i\sigma}(\mathbf{k})$
[$c^{nl}_{i\sigma}(\mathbf{k})$] is the creation [annihilation] operator for
the $l^{\rm th}$ eigenstate in the basis of Wannier functions (see Sec. II.2),
with eigenvalue $E^{l}_{\sigma}(\mathbf{k})$. The summation in Eq. (3) runs
over all Ni(3$d$) and O(2$p$) orbitals; since $n$ and $n^{\prime}$ only
consist of Ni($3d$) orbitals, this corresponds to projecting the Wannier
functions onto Ni($3d$) states, thus capturing the overlap between Ni($3d$)
and O($2p$).

The magnetic moment of a selected nickel site can also be evaluated as

$\displaystyle m_{i}=-\frac{1}{\pi}{\rm
Im}\int\limits_{-\infty}^{E_{F}}d\epsilon$
$\displaystyle\,\mathrm{Tr}_{n}\left[G^{nn}_{ii\uparrow}(\epsilon)-G^{nn}_{ii\downarrow}(\epsilon)\right].$
(4)

Eq. (4) includes in the Ni moments the contribution from the overlap with
oxygen. If we define a different set of Green’s functions by allowing $n$ and
$n^{\prime}$ in Eq. (3) to run over all of the Wannier functions (not just
nickel $3d$), then Eq. (4) corresponds to the total magnetic moment per site,
equivalent to integrating over the density of states near the Fermi level; by
limiting the trace over just the Ni(3$d$) or just O(2$p$) Wannier functions,
we can obtain separate moments on the nickels and oxygens, which may be
compared to those integrated within the PAW spheres Kresse and Joubert (1999),
as discussed in the next section.

###  II.2 Computational approach

We calculated the exchange couplings in Eq. (2) using DFT with the
generalized-gradient approximation exchange-correlation functional of Perdew,
Burke, and Ernzerhof Perdew  _et al._ (1996) and the projector augmented wave
(PAW) methodBlöchl (1994) as implemented in the Vienna ab initio Simulation
Package (VASP) Kresse and Furthmüller (1996a, b). For the rare-earth atoms, we
used PAW potentials corresponding to a $3+$ valence state with $f$-electrons
frozen into the core; for Ni, the 3$p$ semi-core states were included as
valence electrons. We set the energy cutoff of the plane-wave basis to 550 eV
and the energy convergence criteria to 10-8 eV.

Correlation effects were taken into account on the static mean-field level
using the DFT+$U$ method (Anisimov  _et al._ , 1991), where we considered an
effective on-site Coulomb $U$ and Hund’s exchange $J_{H}$ interactions in the
rotationally invariant form Liechtenstein  _et al._ (1995). A value of
$J_{H}=1$ eV can be used for compounds with 3$d$ orbitals Badrtdinov  _et al._
(2016, 2019), while the on-site Coulomb parameter depends on the specific
system. Previous results showed that the best agreement with experiment can be
obtained using a relatively small value of $U$ = 2 eV Hampel and Ederer
(2017); Varignon  _et al._ (2017); Mercy  _et al._ (2017). Therefore, we used
$U$ = 2 eV and $J_{\text{H}}$ = 1 eV for the majority of our DFT+$U$
calculations, and comment on how a larger $U$ ($U$ = 5 eV) may effect the
results.

From the calculated electronic structure, maximally-localized Wannier
functions Marzari and Vanderbilt (1997) were generated using the wannier90
package Pizzi  _et al._ (2020). We constructed Wannier functions for all
occupied and unoccupied Ni 3$d$ and O 2$p$ states. This ensurs that the
obtained Wannier functions for both spin channels are well localized (Fig. 3),
and we can use them as an atomic-like basis for Eqs. (2) and (3). Note, that
this choice is different to Ref. Varignon  _et al._ , 2017, were Wannier
functions were only constructed for the occupied states. Since the $f$
electrons of the RE elements were frozen into core states, we have neglected
their magnetic ordering which only occurs at very low temperatures Fernández-
Díaz  _et al._ (2001); Muñoz  _et al._ (2009).

To examine different magnetic orders we considered two types of cells: the
primitive monoclinic cell containing 20 atoms, and a cell enlarged in the $a$
and $c$ directions containing 80 atoms (necessary for the AFM orderings); an
8$\times$8$\times$6 and 4$\times$6$\times$4 Monkhorst-Pack mesh for Brillouin-
zone integration were used for the 20 and 80 atom cells, respectively. The
number of $k$-points within the first Brillouin zone for the evaluation of the
Green’s function [Eq. (3)] and energy points for integration over a
semicircular contour in the upper half of the complex energy plane [Eqs. (2)
and (4)] were chosen to reach numerical convergence of exchange couplings
within 0.1 meV.

For a systematic study of structural distortion, we used a symmetry-based mode
decomposition analysis Perez-Mato  _et al._ (2010) where we varied the
amplitude of specific modes expressed in terms of irreducible representations
of the parent (cubic perovskite) $Pm\overline{3}m$ structure. The most
significant modes that are present in the RE nickelates include octahedral
rotation modes $M_{3}^{+}$ (in-phase) and $R_{4}^{+}$ (out-of-phase), and the
breathing mode distortion $R_{1}^{+}$ Balachandran and Rondinelli (2013). The
breathing mode, which is responsible for the formation of SB and LB octahedral
(Fig. 1), is the distortion most strongly coupled to the magnetic structure,
so it will be the focus of this study. The main effect of the different RE
atoms is to change the amplitude of the octahedral rotation modes.
Specifically, smaller radii REs correspond to larger rotations. Thus the range
in RE atoms in this study demonstrates the influence of octahedral rotations
on the magnetic ordering. Hampel  _et al._ (2019).

In our calculations, we start from the high-temperature orthorhombic $Pbnm$
structure, relaxed with non-magnetic DFT (with $U=0$) Hampel and Ederer
(2017); in this case, SB and LB octahedra have the same volume. We then
gradually increase the $R_{1}^{+}$ amplitude, calculating how the exchange
couplings and magnetic moments are modified by the distortion. Thus we can
study the evolution of magnetic characteristics with the transition from
$Pbnm$ to the low temperature insulating monoclinic phase $P2_{1}/n$.

Figure 2: (a) T-type antiferromagnetic (AFM) and (b) noncolinear (NCL)
magnetic configurations as a possible ground state magnetic orderings in
$\mathcal{R}$NiO3 with experimental propagation vector
$\mathbf{q}=(\frac{1}{4},\frac{1}{4},\frac{1}{4})$. Yellow and blue spheres
correspond to Ni${}_{\text{SB}}$ and Ni${}_{\text{LB}}$ atoms respectively,
which have smaller and larger magnetic moments corresponding to the length of
the red arrows.

##  III Results

###  III.1 Magnetic moment analysis

We begin with the magnetic moments of PrNiO3, SmNiO3 and LuNiO3, assuming FM
or T-AFM ordering, using the on-site Green’s functions [Eq. (4)] in the
Wannier basis. This procedure will correctly capture the hybridization effect
between nickel and oxygen states and the corresponding distribution of
magnetic moment density Varignon  _et al._ (2017).

Figure 3:  Maximally localized Wannier functions of the long-bond (LB) and
short-bond (SB) nickel atoms. (a) corresponds to $d_{x^{2}-y^{2}}$ and (b) to
$d_{z^{2}}$ orbitals for nearest neighbor LB – SB – LB nickel atoms,
calculated for the T-AFM PrNiO3 system at the experimental breathing mode
amplitude $R_{1}^{+}$ = 0.045 Å. Only spin majority orbitals at each nickel
atoms are represented for clarity, spin moments are shown in (c).

Nickel magnetic moments versus $R_{1}^{+}$ amplitude are shown in Fig.
4(a)-(c). Dashed vertical lines denote the experimental breathing-mode
amplitudes for PrNiO3 (0.045 ÅMedarde  _et al._ (2008)) and LuNiO3 (0.075
ÅAlonso  _et al._ (2001)), and the theoretically predicted amplitude for
SmNiO3 (0.060 ÅHampel  _et al._ (2019)), where no experimental data on the low
temperature phase is available. In the case of the high-temperature
orthorhombic phase ($R_{1}^{+}$ = 0 Å), in which all Ni octahedra have the
same volume, FM order leads to magnetic moments of $\simeq 0.9\,\mu_{B}$ per
nickel atom for all materials. Within the T-AFM order the $Pbnm$ symmetry is
broken, causing a small variation of magnetic moments for the inequivalent Ni
sites. The difference between the SB and LB moments strongly increases with
larger $R_{1}^{+}$ breathing mode amplitudes. In the FM case, the SB Ni
magnetic moment smoothly decreases; for T-AFM, the moment quickly vanishes at
a modest value of $R_{1}^{+}$. The amplitude where the SB moment goes to zero
depends on the rare-earth element, specifically it is larger for LuNiO3 than
PrNiO3 and SmNiO3. In all cases, the $R_{1}^{+}$ amplitudes at which the SB
moment vanish are much smaller than the experimentally observed amplitudes for
these compounds. Therefore, DFT+$U$ predicts $m_{\rm SB}=0$ for T-AFM order
for all experimental structures. In contrast, $m_{\text{LB}}$ increases with
$R_{1}^{+}$ and converges to $\simeq 1.3\,\mu_{B}$ for both FM and T-AFM
order.

The Wannier-based analysis allows us to analyze the orbital contributions to
total nickel magnetic moments, and we find that for both FM and T-AFM order,
LB and SB moments originate almost entirely from Ni $e_{g}$ states. Since the
corresponding Ni $t_{2g}$ states are fully occupied, they are not magnetically
active.

Figure 4: Evolution of the SB and LB magnetic moments of (a) PrNiO3, (b)
SmNiO3, and (c) LuNiO3 as a function of the $R_{1}^{+}$ amplitude, calculated
for different magnetic states using $U$ = 2 eV. Dashed vertical lines
correspond to experimentally observed breathing mode amplitudes from Ref. 45
for PrNiO3 and Ref. 46 for LuNiO3, and the theoretically predicted amplitude
for SmNiO3 from Ref. 44. The $R_{1}^{+}$ = 0 structure corresponds to the high
temperature orthorhombic $Pbnm$ structure. Additional dotted curves for PrNiO3
denoted by crosses represent results obtained with $U$ = 5 eV.

These obtained results are consistent with the picture that would emerge from
an analysis of the magnetization within the PAW spheres around Ni and O. The
nickel moment from the Wannier analysis (and its dependence on $R_{1}^{+}$) is
in good agreement ($\sim 0.01$ $\mu_{B}$) with the moments in Ni PAW spheres.
This indicates that the magnetism is localized around the atoms, instead of in
the interstitial regions Badrtdinov  _et al._ (2018). Similarly to the PAW
estimation, the Wannier based oxygen magnetic moments for all magnetic
configurations are very small (less than 0.05 $\mu_{B}$). This suggests that,
while there is significant Ni-O hybridization, the moment resides only at the
Ni sites.

The magnitude of the magnetic moments depend on the choice of $U$ parameter of
DFT+$U$. Increasing $U$ leads to a stronger localization of the $3d$
electrons, and consequently gives comparably larger total magnetic moments of
nickel atoms. The reason for this can be seen in the density of states, e.g.,
for PrNiO3 (Fig. 5), where a larger $U$ value leads to stronger splitting
between spin channels of LB nickel $e_{g}$ states, enhancing their magnetic
moment. To demonstrate this effect, we show magnetic moments calculated using
$U$ = 5 eV for PrNiO3 [see crosses in Fig. 4(a)]. The resulting magnetic
moments for $Pbnm$ are larger than those obtained with $U$ = 2 eV. For
instance, the large $R_{1}^{+}$ amplitude LB moments (for both magnetic
orders) converge to $\simeq 1.3\,\mu_{B}$ for $U$ = 2 eV versus $\simeq
1.5\,\mu_{B}$ for $U$ = 5 eV. Furthermore, $m_{\text{SB}}$ for T-AFM order
becomes zero at comparably larger value $R_{1}^{+}$ $\simeq$ 0.05 Å compared
to $\leq$ 0.007 Å for $U$ = 2 eV, as was found in previous DFT+$U$ studies
Hampel and Ederer (2017).

Figure 5: Partial densities of states (DOS) for PrNiO3 with the experimental
$R_{1}^{+}$ = 0.045 Å, calculated with T-AFM order for (a) $U$ = 2 eV and (b)
$U$ = 5 eV. Here, only $t_{2g}$ and $e_{g}$ states contributions of a single
LB and SB nickel atom, as well as 2$p$ states of the surrounding six oxygen
atoms in NiO6 octahedra are shown. The different occupation of the $e_{g}$
spin channels of the LB nickel site carries the magnetic moment, while for SB
both spin channels are equally occupied, leading to a zero moment.

###  III.2 Exchange couplings

We now perform a quantitative analysis of the various magnetic interactions
present in the $\mathcal{R}$NiO3 materials with respect to $R_{1}^{+}$
amplitude. In Fig. 6 we show a schematic of the exchange interactions included
in our study. The LB and SB Ni sublattices are depicted as blue and yellow
spheres, and the RE and oxygens are removed for clarity. Focusing on a LB
site, there are two inequivalent nearest-neighbor interactions with SB sites
in the $ab$ plane ($J_{1}$). Also in the $ab$ plane, we consider next-nearest-
neighbor interactions with other LB sites $J_{2}$ and $J_{3}$, as well as
next-next-nearest neighbor interactions $J_{4}$. In the $c$ direction, we have
nearest neighbor interactions between LB and SB sites ($J_{5}$), next-nearest
neighbor ($J_{7}$ and $J_{6}$) interactions.

Figure 6:  Magnetic model of $\mathcal{R}$NiO3 system with main exchange
interactions plotted in (a) side and (b) top view. Blue atoms are long bond
(LB) Ni, yellow are short bond (SB) (all other atoms removed). Arrows
represent magnetic moments of nickel atoms in the T-AFM order. Corresponding
interactions between SB nickels are also included, but not shown in the figure
for clarity.

####  III.2.1 Ferromagnetic interactions $J_{1}$ and $J_{5}$

Figure 7: Nearest neighbor (a) $J_{1}^{a/b}$ and (b) $J_{5}$ exchange
couplings (see Fig. 6) between nickel atoms in $\mathcal{R}$NiO3 as a function
of breathing mode distortion $R_{1}^{+}$ amplitude. The reference magnetic
configuration is FM. The additional dotted curves for PrNiO3 denoted by
crosses represent results obtained with $U$ = 5 eV.

The calculated magnetic interactions via Eq. (2) depend on the reference
magnetic configuration used in the underlying DFT calculation Zhu  _et al._
(2020), since the electronic structure and the resulting single-particle
Green’s functions will be different for FM and T-AFM order. Since the T-AFM
order leads to a suppression of the SB magnetic moment (Fig. 4) the nearest-
neighbor exchange interactions $J_{1}$ and $J_{5}$ are suppressed, as they
couple two different types of nickel sites in the $ab$ plane and along the $c$
direction (Fig. 6). Therefore, for these interactions we focus on the FM state
as the reference magnetic order. Resulting values as a function of $R_{1}^{+}$
are given in Fig. 7.

In general, short range interactions originate from a combination of different
mechanisms Badrtdinov  _et al._ (2016, 2019). Previous studies reported that
the charge fluctuation between SB Ni and LB Ni sites with unequal spins will
lead to ferromagnetic double exchange $J^{\text{DE}}<0$ process, which
competes with the conventional antiferromagnetic superexchange
$J^{\text{SE}}>0$ Lu  _et al._ (2018). The resulting coupling
$J=J^{\text{SE}}+J^{\text{DE}}$ can have either sign depending on which
mechanism will prevail. We show in Fig. 7 strong FM interaction (i.e.,
negative values) for both $J_{1}$ and $J_{5}$, implying
$|J^{\text{SE}}|<|J^{\text{DE}}|$ for both cases. Double exchange between the
Ni atoms requires inequivalent charges on the sites Zener (1951); Keshavarz
_et al._ (2017). The strength of double exchange, and thus the FM couplings,
would be expected to increasing with $R_{1}^{+}$ Lu  _et al._ (2018). However,
in our calculations the magnitude of $J_{1}$ and $J_{5}$ couplings gradually
_decrease_ with increasing $R_{1}^{+}$ amplitudes. This suggests that other
factors (e.g., charge transfer to oxygen states) are important in the
resulting exchange interactions.

As a result of the presence of $M_{3}^{+}$ and $R_{4}^{+}$ rotation modes, the
nearest neighbor exchange $J_{1}$ splits into two groups $J^{a}_{1}$ and
$J^{b}_{1}$ for finite $R_{1}^{+}$ (Fig. 6). The value of this splitting
depends on the rare-earth element $\mathcal{R}$ (c.f., solid and dashed curves
in Fig. 7). In particular, a stronger splitting $J^{a}_{1}<J^{b}_{1}$ can be
found for LuNiO3 compared to PrNiO3, where $J^{a}_{1}\simeq J^{b}_{1}$. Hence,
this magnitude of the splitting is connected to the magnitude of the
octahedral rotation modes $M_{3}^{+}$ and $R_{4}^{+}$ Hampel and Ederer
(2017). As mentioned in Sec. II.2, RE elements with smaller effective radii
will cause stronger $M_{3}^{+}$ and $R_{4}^{+}$ rotations of NiO6 octahedra,
modifying Ni – O – Ni angles, and resulting in a larger difference between
$J^{a}_{1}$ and $J^{b}_{1}$.

Furthermore, higher $U$ parameters will lead to a stronger localization of $d$
electrons, which influences not only the magnetic moments (Fig. 4), but also
the corresponding exchange interactions between them. Our calculations on
PrNiO3 with $U$ = 5 eV show an enhancement of the nearest neighbor
ferromagnetic $J_{1}$ and $J_{5}$ couplings (Fig. 7). This is directly related
to the fact that using a larger on-site Coulomb interaction parameter will
decrease the AFM superexchange term, which is inversely proportional to $U$.
Double exchange is proportional to the transfer integral between Ni site, and
does not depend directly on $U$ Keshavarz  _et al._ (2017).

####  III.2.2  Antiferromagnetic interactions $J_{4}$ and $J_{7}$

Figure 8:  Next nearest neighbor $J_{4}^{a/b}$ magnetic interactions between
nickel atoms in $\mathcal{R}$NiO3 as a function of breathing mode distortion
$R_{1}^{+}$ and reference magnetic order: (a) interactions between LB atoms in
T-AFM order, while (b) and (c) correspond to interactions between LB and SB in
FM order, respectively. Notation of magnetic couplings corresponds to Fig. 6.
Figure 9:  Next nearest neighbor $J_{7}$ magnetic interactions between nickel
atoms in $\mathcal{R}$NiO3 as a function of breathing mode distortion
$R_{1}^{+}$ and reference magnetic order: (a) interactions between LB atoms in
T-AFM order, while (b) and (c) correspond to interactions between LB and SB in
FM order, respectively. Additional dotted curve for PrNiO3 denoted by crosses
represents results obtained with $U$ = 5 eV. Notation of magnetic couplings
corresponds to Fig. 6.

Next-nearest-neighbor exchange interactions $J_{4}$ and $J_{7}$ connect the
same type (SB or LB) of nickel atom (Fig. 6). Similar to the nearest-neighbor
$J_{1}$ interaction, the next nearest neighbor $J_{4}$ also splits into two
groups $J^{a}_{4}$ and $J^{b}_{4}$. The value of this splitting depends on the
RE element via the octahedral rotation modes $M_{3}^{+}$ and $R_{4}^{+}$.
These interactions are antiferromagnetic, i.e., positive (Figs. 8 and 9), and
the dependence on the breathing mode distortion is different depending on the
reference magnetic order.

For the case of FM order, the next-nearest neighbor exchange interactions
between LB nickel atoms gradually increase, while for the SB nickels they
gradually decrease to zero; this is the same trend as the magnetic moments,
see Fig. 4. In contrast, for the T-AFM order, the SB magnetic moments are
suppressed and thus the magnetic interactions between LB nickel atoms jump
abruptly to a large value with $R_{1}^{+}>0$. Further increasing $R_{1}^{+}$
results in a small reduction in the magnitude of the exchange coupling. The
behavior in Fig. 8(a) hints at a strong coupling between $R_{1}^{+}$
distortion and T-AFM order. Overall, exchange couplings between LB nickel
sites at both reference magnetic configurations have the same order of
magnitude. The values of $J_{4}$ and $J_{7}$ at $R_{1}^{+}$ = 0 are about an
order of magnitude smaller than $J_{1}$/$J_{5}$, and thus the T-AFM magnetic
configuration is fragile. Since the main mechanism of the next nearest
interactions is superexchange, all these couplings weaken for a larger $U$
value, which is shown for $J_{7}$ of PrNiO3 in Fig. 9.

####  III.2.3 Summary and remaining exchange couplings

Table 1: The main exchange couplings (in meV) for PrNiO3, SmNiO3 and LuNiO3, calculated ($U$=2 eV) using experimental breathing mode amplitudes $R_{1}^{+}$ for PrNiO3 and LuNiO3, and the theoretical one for SmNiO3. Notation of magnetic couplings corresponds to Fig. 6.  | PrNiO3 | SmNiO3 | LuNiO3
---|---|---|---
$R_{1}^{+}$ | 0.045 Å Medarde  _et al._ (2008) | 0.060 Å Hampel  _et al._ (2019) | 0.075 Å (Alonso  _et al._ , 2001)
$J_{1a}/J_{1b}$ | -103.9/-102.6 | -89.4/-81.4 | -83.6/-59.8
$J_{5}$ | -96.7 | -81.4 | -60.7
$J_{4a}/J_{4b}$ | 38.7/37.5 | 33.0/28.3 | 29.8/17.9
$J_{7}$ | 41.3 | 26.9 | 7.7
$J_{2}$ | 3.5 | 7.5 | 10.0
$J_{3}$ | 9.9 | 11.6 | 16.0
$J_{6a}/J_{6b}$ | 10.9/7.2 | 7.8/7.7 | 6.2/8.5

The values for the exchange interactions discussed above, calculated at the
experimental $R_{1}^{+}$ amplitudes Medarde  _et al._ (2008); Alonso  _et al._
(2001) (except for SmNiO3 where we use the theoretical value from Ref. Hampel
_et al._ , 2019), are given in Table 1. Also included are the remaining
couplings shown in Fig. 6: $J_{2}$, $J_{3}$, and $J_{6}$; we see that these
couplings are small compared to $J_{1}/J_{5}$ and $J_{4}/J_{7}$. This is
because of the geometry of the active $e_{g}$ orbitals, which overlap with
each other via oxygen atoms along the Cartesian axes (Fig. 6). Our orbital
resolved analysis of exchange interactions (Eq. 2) demonstrates that the main
contribution to in-plane $J_{1}$ and $J_{4}$ comes from $d_{x^{2}-y^{2}}$
orbitals, while out-of-plane $J_{5}$ and $J_{7}$ are originated from
$d_{z^{2}}$ states. Furthermore, it can be seen that the exchange couplings
significantly increase going from Lu to Pr, e.g. when octahedral rotations are
reduced.

###  III.3 Anisotropy and non-collinear order

To shed light on the long-standing question of collinear versus noncollinear
(NCL) order in rare-earth nickelates (Fig. 2), we also performed an analysis
of the magnetic anisotropy and NCL magnetic arrangements within DFT+$U$
calculations taking spin-orbit coupling (SOC) into account. In the case of the
proposed Scagnoli  _et al._ (2006); Lu  _et al._ (2018) NCL magnetic order,
nearest-neighbor magnetic moments are orthogonal to each other [Fig. 2(b)].
Since $\mathbf{m}_{\rm LB}\cdot\mathbf{m}_{\rm SB}=0$, the contributions of
the strong isotropic ferromagnetic $J_{1}/J_{5}$ couplings are eliminated, and
magnetic frustration between FM nearest-neighbors and AFM next-nearest-
neighbors couplings does not occur.

We perform calculations including SOC on the $Pbnm$ crystal structure of
PrNiO3, SmNiO3, and LuNiO3 with breathing mode amplitudes given in Table 1. We
find that the SB moments for all three materials vanishes in these
calculations, analogously to the collinear calculations in Fig. 4. The LB
moment is always nonzero and equals 1.19 $\mu_{B}$, 1.24 $\mu_{B}$ and 1.28
$\mu_{B}$ for $\mathcal{R}$ = Pr, Sm and Lu, respectively. Thus, the NCL
configuration we find is equivalent to the collinear T-AFM order, with the LB
magnetic moments free to align along the direction which minimizes the
anisotropy energy in the systems; we can find this direction by calculating
the anisotropic energy map Abdeldaim  _et al._ (2019).

Figure 10: Anisotropic energy map for PrNiO3 system. The arrow indicates the
energy minimum, which is close to the direction of apical oxygen position of
MoO6 octahedra, depicted in (b) together with Wannier function of
$x^{2}-y^{2}$ symmetry, the valence electron in which state induces orbital
magnetic moment.

In Fig. 10(a), we show the interpolated anisotropic energy map for PrNiO3. To
generate this plot, DFT+$U$+SOC calculations were performed using magnetic
moments constrained along various (uniformly distributed) directions, but
keeping their mutual orientation consistent with T-AFM order. We find an easy
axis along a direction close to the position of the apical oxygen atom in the
NiO6 octahedra [green arrow in Fig. 10(a)]. The anisotropy in the energy is
very small ($\sim$ 0.03 meV), which is expected for $3d$ systems with weak SOC
Badrtdinov  _et al._ (2016, 2019). According to conventional theories of
magnetocrystalline anisotropy (i.e., Bruno’s model Bruno (1989); Bru (2006)),
the orbital moment is maximized when the magnetization points along the easy
axis. This is true for $\mathcal{R}$NiO3 systems, showing the maximum value of
orbital moment $M_{L}$ = 0.13 $\mu_{B}$ aligned with spin magnetic moment 1.19
$\mu_{B}$ along the direction of the easy axis.

The orbital moment is induced by $e_{g}$ electrons, which can be seen on the
level of density matrix $\hat{N}$:

$\displaystyle M_{L}^{z}={\rm Tr}_{m}(\hat{N}\cdot\hat{L}_{z}).$

In particular, the valence electron occupying the $d_{z^{2}}$ will not produce
any orbital moment ($\hat{L}_{z}d_{z^{2}}=0$). On the other hand, due to a
non-zero $\braket{d_{xy}}{\hat{L}_{z}}{d_{x^{2}-y^{2}}}=2i$, and the fact that
the calculated density matrix contains imaginary components between $d_{xy}$
and $d_{x^{2}-y^{2}}$ states, an orbital moment is induced along the direction
of the apical oxygen in the NiO6 octahedra ($z$-axis) by the electron in the
$d_{x^{2}-y^{2}}$ state. This leads to the formation of the easy axis
anisotropy found by in our DFT+$U$+SOC calculations (Fig. 10).

##  IV Discussion

The goal of the calculations in Sec. III was to elucidate two key issues
regarding the magnetic structure in RE nickelates: (a) why DFT+$U$ predicts a
FM structure, while an AFM structure with propagation vector
$\mathbf{q}=(\frac{1}{4},\frac{1}{4},\frac{1}{4})$ is observed experimentally;
and (b) what is the relation between the T-AFM magnetic structure determined
from neutron diffraction data, and the NCL magnetic structure determined from
RIXS.

We find that the key physical mechanism at play in these systems is the
competition between strong nearest-neighbor ferromagnetic $J_{1}/J_{5}$ and
the much weaker next-nearest-neighbor antiferromagnetic $J_{4}/J_{7}$. This
competition causes magnetic frustration, and works to destabilize the long-
range magnetic order. As discussed before, the $R_{1}^{+}$ distortion will
result in different magnetic moments on the LB and SB Ni atoms, with the SB
magnetic moment suppressed (Fig. 4). Therefore, the energy gain from forming
ferromagnetic order will be directly proportional to the SB moment. Hence, the
reason that DFT+$U$ favors FM order Hampel and Ederer (2017); Varignon  _et
al._ (2017) is the overestimation of short bond magnetic moments, leading to a
larger contribution of nearest-neighbor FM couplings, which makes the FM
solution energetically favorable. For a sufficiently small SB moment the FM
$J_{1}$ and $J_{5}$ interactions are suppressed, and the next-nearest neighbor
AFM couplings $J_{4}$ and $J_{7}$ will dominate; the result is the collinear
T-AFM order proposed from neutron diffraction experiments García-Muñoz  _et
al._ (1994); Alonso  _et al._ (1999); Muñoz  _et al._ (2009); Gawryluk  _et
al._ (2019) (Fig. 2). As we saw in Fig. 4, the SB moments go to zero rapidly
with increasing $R_{1}^{+}$ amplitude in the T-AFM ordering to avoid the
energy penalty associated with a finite FM $J_{1}/J_{5}$ coupling.

The reason for the theory-experiment discrepancy, could be due to an
overestimation of exchange splitting in DFT+$U$, or problems in the refinement
of experimental data. For a future analysis it would be of high interest to
incorporate self-energy corrections to the Green’s function exchange coupling
analysis from higher order theory, i.e, GW or dynamical mean-field theory, to
rule out problems due to dynamic correlation effects Kvashnin  _et al._
(2015).

Another approach to avoiding frustration is the proposed magnetic structure
from RIXS experiments Scagnoli  _et al._ (2006); Lu  _et al._ (2018). The NCL
ordering where the nearest-neighbor moments point in orthogonal directions,
has no energy penalty for AFM versus FM ordering, and thus does not require
the strong reduction of the SB magnetic moment needed for stabilization of
collinear antiferromagnetic order. However, as discussed in Sec. III.3, we
cannot stabilize this ordering in our calculations, as the SB moment vanishes
at the experimental $R_{1}^{+}$ amplitude, even if a SOC is included (thus a
NCL order is allowed). The origin of the SB magnetic suppression can be
related to the symmetry of the magnetic configurations. The SB magnetic moment
is surrounded by four LB nickel atoms within the $ab$ plane (Fig. 2), and in
the case of FM order, this will enhance local magnetic polarization on the SB
nickel atom, increasing its magnetic moment. But since in T-AFM (and also the
NCL order) half of the neighboring LB moments are aligned parallel, and the
other half anti-parallel, the SB moment will be reduced. The valence electron
in the magnetically active $d_{x^{2}-y^{2}}$ states induces an orbital moment
$M_{L}\simeq$ 0.13 $\mu_{B}$ pointed at the apical oxygen position of the NiO6
octahedra. This leads to the formation of easy-axis anisotropy, giving
preferable direction for LB magnetic moments alignment.

The analysis above is for the experimental $R_{1}^{+}$ amplitudes (theoretical
for SmNiO3); however, we see that these results are robust over a wide range
of $R_{1}^{+}$ amplitudes. Increasing the breathing distortion reduces FM
nearest neighbor $J_{1}$/ $J_{5}$ interactions and magnetic moments of SB
nickel atoms, thus it can strengthen long range T-AFM ordering in these
systems. The same behavior was found in previous calculations Hampel and
Ederer (2017) by directly comparing relaxed $R_{1}^{+}$ amplitudes in DFT+$U$,
which revealed $R_{1}^{+}$ (FM) $<$ $R_{1}^{+}$(T-AFM).

The experimental estimates, used to fit magnon branches in RIXS experiment for
thin film NdNiO3 samples, revealed that the in $ab$-plane next-nearest-
neighbor interaction $J_{4}$ is stronger that the nearest neighbor $J_{2}$
coupling $J_{4}\simeq 2J_{2}$ Lu  _et al._ (2018). In our calculations we also
observe $J_{4}\gg J_{2}$ (Table 1). Furthermore, we found strong spacial
anisotropy depending on RE element, which splits interactions in the $ab$
plane into two groups: nearest-neighbor $J_{1a}/J_{1b}$, next-nearest-neighbor
$J_{2}/J_{3}$ and next-next-nearest neighbor $J_{4a}/J_{4b}$ (Fig. 6).
Moreover, we see a strong difference between $ab$ plane ($J_{1}/J_{4}$) and
out-of-plane $c$-direction ($J_{5}/J_{7}$) exchange couplings depending on RE
element.

Finally, the results provide a clear picture of the trends of the magnetic
order across the nickelate series. While the larger octhadral distortions
($R_{4}^{+},M_{3}^{+}$) present for smaller RE ions make the system more
susceptible to the insulating breathing-mode-distorted state Peil  _et al._
(2019); Mercy  _et al._ (2017), we find that they also reduce the exchange
couplings, making the formation of the magnetically ordered state less
favorable. This is in agreement with the experimentally-observed reduction in
stability of magnetic order in the low temperature insulating phase across the
series towards LuNiO3 Hampel and Ederer (2017); Catalano  _et al._ (2018);
Gawryluk  _et al._ (2019).

##  V Conclusions

Based on DFT+$U$ calculations, Wannier orbitals and Green’s function
techniques, we constructed microscopic magnetic models for selected rare-earth
nickelate systems (PrNiO3, SmNiO3, and LuNiO3) and studied the evolution of
their parameters under the $R_{1}^{+}$ breathing mode distortion and reference
magnetic configuration. Our theoretical results help to clarify several issues
related to these systems. The Wannier function based analysis demonstrates
that the magnetic moments reside on nickel atoms, and they are induced by
active Ni $e_{g}$ states. The main exchange interactions are strong
ferromagnetic interactions arising from double exchange between the nearest
neighbor nickels, and somewhat weaker antiferromagnetic interactions between
the next-nearest neighbor nickels induced by a kinetic superexchange
mechanism. The competition between these, together with the $R_{1}^{+}$
distortion, determines the magnetic order.

The stronger nearest-neighbor interactions are the reason that DFT+$U$
calculations find the ferromagnetic state to be lowest in energy. The
stabilization of collinear T-AFM order, which is compatible with the
experimental symmetry, is possible due to suppression of short bond magnetic
moments with the onset of the breathing mode distortion, eliminating the
contribution of the ferromagnetic interactions. We find that even if spin-
orbit coupling is included and a non-collinear ordering is allowed, the SB Ni
moments are still suppressed, and T-AFM order is restored with an easy axis
magnetization pointed towards the apical oxygen. Thus our work determines the
reason for the general disagreement between the DFT-$U$ predicted magnetic
order (FM) and experiment (AFM), as well as provides a clear explanation how
the $\mathbf{q}=(\frac{1}{4},\frac{1}{4},\frac{1}{4})$ magnetic order emerges
in rare-earth nickelates, and how the trends evolve across the series.

###### Acknowledgements.

D.I.B. is grateful to S. A. Nikolaev (Tokyo Institute of Technology) for
helpful discussions. The calculations have been performed using the facilities
of the Flatiron Institute. The Flatiron Institute is a division of the Simons
Foundation. C.E.D. acknowledge support from the National Science Foundation
under Grant No. DMR-1918455.

## References

  * Hampel and Ederer (2017) A. Hampel and C. Ederer, Phys. Rev. B 96, 165130 (2017).
  * Varignon  _et al._ (2017) J. Varignon, M. N. Grisolia, J. Íñiguez, A. Barthélémy, and M. Bibes, npj Quantum Materials 2, 21 (2017).
  * Catalano  _et al._ (2018) S. Catalano, M. Gibert, J. Fowlie, J. Íñiguez, J.-M. Triscone, and J. Kreisel, Rep. Prog. Phys. 81, 046501 (2018).
  * Lu  _et al._ (2017) Y. Lu, Z. Zhong, M. W. Haverkort, and P. Hansmann, Phys. Rev. B 95, 195117 (2017).
  * Obradors  _et al._ (1993) X. Obradors, L. M. Paulius, M. B. Maple, J. B. Torrance, A. I. Nazzal, J. Fontcuberta, and X. Granados, Phys. Rev. B 47, 12353 (1993).
  * Catalano  _et al._ (2014) S. Catalano, M. Gibert, V. Bisogni, O. E. Peil, F. He, R. Sutarto, M. Viret, P. Zubko, R. Scherwitzl, A. Georges, G. A. Sawatzky, T. Schmitt, and J.-M. Triscone, APL Materials 2, 116110 (2014).
  * Boris  _et al._ (2011) A. V. Boris, Y. Matiks, E. Benckiser, A. Frano, P. Popovich, V. Hinkov, P. Wochner, M. Castro-Colin, E. Detemple, V. K. Malik, C. Bernhard, T. Prokscha, A. Suter, Z. Salman, E. Morenzoni, G. Cristiani, H.-U. Habermeier, and B. Keimer, Science 332, 937 (2011).
  * Y  _et al._ (2015) C. Y, B. L., Y. Sun, H. Yang, and X. Zhang, Electrochimica Acta 174, 41 (2015).
  * Shi  _et al._ (2014) J. Shi, Y. Zhou, and S. Ramanathan, Nature Communications 5, 4860 (2014).
  * Cheong and Mostovoy (2007) S.-W. Cheong and M. Mostovoy, Nature Materials 6, 13 (2007).
  * Alonso  _et al._ (1999) J. A. Alonso, J. L. García-Muñoz, M. T. Fernández-Díaz, M. A. G. Aranda, M. J. Martínez-Lope, and M. T. Casais, Phys. Rev. Lett. 82, 3871 (1999).
  * Mazin  _et al._ (2007) I. I. Mazin, D. I. Khomskii, R. Lengsdorf, J. A. Alonso, W. G. Marshall, R. M. Ibberson, A. Podlesnyak, M. J. Martínez-Lope, and M. M. Abd-Elmeguid, Phys. Rev. Lett. 98, 176406 (2007).
  * Johnston  _et al._ (2014) S. Johnston, A. Mukherjee, I. Elfimov, M. Berciu, and G. A. Sawatzky, Phys. Rev. Lett. 112, 106404 (2014).
  * Green  _et al._ (2016) R. J. Green, M. W. Haverkort, and G. A. Sawatzky, Phys. Rev. B 94, 195127 (2016).
  * Park  _et al._ (2012) H. Park, A. J. Millis, and C. A. Marianetti, Phys. Rev. Lett. 109, 156402 (2012).
  * Lau and Millis (2013) B. Lau and A. J. Millis, Phys. Rev. Lett. 110, 126404 (2013).
  * Subedi  _et al._ (2015) A. Subedi, O. E. Peil, and A. Georges, Phys. Rev. B 91, 075128 (2015).
  * Haule and Pascut (2017) K. Haule and L. Pascut, G, Scientific Reports 7, 10375 (2017).
  * Peil  _et al._ (2019) O. E. Peil, A. Hampel, C. Ederer, and A. Georges, Phys. Rev. B 99, 245127 (2019).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * García-Muñoz  _et al._ (1994) J. L. García-Muñoz, J. Rodríguez-Carvajal, and P. Lacorre, Phys. Rev. B 50, 978 (1994).
  * Gawryluk  _et al._ (2019) D. J. Gawryluk, Y. M. Klein, T. Shang, D. Sheptyakov, L. Keller, N. Casati, P. Lacorre, M. T. Fernández-Díaz, J. Rodríguez-Carvajal, and M. Medarde, Phys. Rev. B 100, 205137 (2019).
  * Muñoz  _et al._ (2009) A. Muñoz, J. A. Alonso, M. J. Martínez-Lope, and M. T. Fernández-Díaz, J. Solid State Chem. 182, 1982 (2009).
  * Fernández-Díaz  _et al._ (2001) M. T. Fernández-Díaz, J. A. Alonso, M. J. Martínez-Lope, M. T. Casais, and J. L. García-Muñoz, Phys. Rev. B 64, 144417 (2001).
  * Scagnoli  _et al._ (2006) V. Scagnoli, U. Staub, A. M. Mulders, M. Janousch, G. I. Meijer, G. Hammerl, J. M. Tonnerre, and N. Stojic, Phys. Rev. B 73, 100409 (2006).
  * Lu  _et al._ (2018) Y. Lu, D. Betto, K. Fürsich, H. Suzuki, H.-H. Kim, G. Cristiani, G. Logvenov, N. B. Brookes, E. Benckiser, M. W. Haverkort, G. Khaliullin, M. Le Tacon, M. Minola, and B. Keimer, Phys. Rev. X 8, 031014 (2018).
  * Liechtenstein  _et al._ (1987) A. Liechtenstein, M. Katsnelson, V. Antropov, and V. Gubanov, Journal of Magnetism and Magnetic Materials 67, 65 (1987).
  * Mazurenko and Anisimov (2005) V. V. Mazurenko and V. I. Anisimov, Phys. Rev. B 71, 184434 (2005).
  * Momma and Izumi (2011) K. Momma and F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011).
  * Note (1) For simplicity, we use spin quantum number $S=\frac{1}{2}$ in Eq. (2) during the whole range of breathing distortions. It was revealed that very large breathing distortion amplitudes are needed to achieve spins $0$ and $1$ for SB and LB magnetic moments Johnston  _et al._ (2014); Green  _et al._ (2016), which is much higher than our considered cases.
  * Perdew  _et al._ (1996) J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Computa. Mater. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
  * Anisimov  _et al._ (1991) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B 44, 943 (1991).
  * Liechtenstein  _et al._ (1995) A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467 (1995).
  * Badrtdinov  _et al._ (2016) D. I. Badrtdinov, O. S. Volkova, A. A. Tsirlin, I. V. Solovyev, A. N. Vasiliev, and V. V. Mazurenko, Phys. Rev. B 94, 054435 (2016).
  * Badrtdinov  _et al._ (2019) D. I. Badrtdinov, V. V. Mazurenko, and A. A. Tsirlin, Phys. Rev. B 100, 214401 (2019).
  * Mercy  _et al._ (2017) A. Mercy, J. Bieder, J. Íñiguez, and P. Ghosez, Nat. Commun. 8, 1677 (2017).
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  * Pizzi  _et al._ (2020) G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, Journal of Physics: Condensed Matter 32, 165902 (2020).
  * Perez-Mato  _et al._ (2010) J. M. Perez-Mato, D. Orobengoa, and M. I. Aroyo, Acta Crystallographica Section A 66, 558 (2010).
  * Balachandran and Rondinelli (2013) P. V. Balachandran and J. M. Rondinelli, Physical Review B 88, 054101 (2013).
  * Hampel  _et al._ (2019) A. Hampel, P. Liu, C. Franchini, and C. Ederer, npj Quantum Materials 4, 5 (2019).
  * Medarde  _et al._ (2008) M. Medarde, M. T. Fernández-Díaz, and P. Lacorre, Phys. Rev. B 78, 212101 (2008).
  * Alonso  _et al._ (2001) J. A. Alonso, M. J. Martínez-Lope, M. T. Casais, J. L. García-Muñoz, M. T. Fernández-Díaz, and M. A. G. Aranda, Phys. Rev. B 64, 094102 (2001).
  * Badrtdinov  _et al._ (2018) D. I. Badrtdinov, S. A. Nikolaev, A. N. Rudenko, M. I. Katsnelson, and V. V. Mazurenko, Phys. Rev. B 98, 184425 (2018).
  * Zhu  _et al._ (2020) X. Zhu, A. Edström, and C. Ederer, Phys. Rev. B 101, 064401 (2020).
  * Zener (1951) C. Zener, Phys. Rev. 82, 403 (1951).
  * Keshavarz  _et al._ (2017) S. Keshavarz, Y. O. Kvashnin, D. C. M. Rodrigues, M. Pereiro, I. Di Marco, C. Autieri, L. Nordström, I. V. Solovyev, B. Sanyal, and O. Eriksson, Phys. Rev. B 95, 115120 (2017).
  * Abdeldaim  _et al._ (2019) A. H. Abdeldaim, D. I. Badrtdinov, A. S. Gibbs, P. Manuel, H. C. Walker, M. D. Le, C. H. Wu, D. Wardecki, S.-G. Eriksson, Y. Kvashnin, A. A. Tsirlin, and G. J. Nilsen, Phys. Rev. B 100, 214427 (2019).
  * Bruno (1989) P. Bruno, Phys. Rev. B 39, 865 (1989).
  * Bru (2006) _Magnetism. From fundamentals to nanoscale dynamics_ (Springer, Berlin Heidelberg, 2006).
  * Kvashnin  _et al._ (2015) Y. O. Kvashnin, O. Grånäs, I. Di Marco, M. I. Katsnelson, A. I. Lichtenstein, and O. Eriksson, Phys. Rev. B 91, 125133 (2015).