# Large dynamical magnetic effective charges and anti-magnetoelectricity from
spin and orbital origin in multiferroic BiCoO3

Maxime Braun  Physique Théorique des Matériaux, QMAT, CESAM, Université de
Liège, B-4000 Sart-Tilman, Belgium  Univ. Lille, CNRS, Centrale Lille, ENSCL,
Univ. Artois, UMR 8181-UCCS-Unité de Catalyse et Chimie du Solide, F-59000
Lille, France  Bogdan Guster  Physique Théorique des Matériaux, QMAT, CESAM,
Université de Liège, B-4000 Sart-Tilman, Belgium  Andrea Urru  Department of
Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854, USA
Houria Kabbour  Univ. Lille, CNRS, Centrale Lille, ENSCL, Univ. Artois, UMR
8181-UCCS-Unité de Catalyse et Chimie du Solide, F-59000 Lille, France  Eric
Bousquet  Physique Théorique des Matériaux, QMAT, CESAM, Université de Liège,
B-4000 Sart-Tilman, Belgium

###### Abstract

Using first-principles calculations, we explore the magnetoelectric properties
of the room-temperature multiferroic crystal BiCoO3. We use both applied
magnetic field and finite-difference techniques to show that BiCoO3 is anti-
magnetoelectric at the linear level. The calculation of the dynamical
effective charges reveals that the total magnetoelectric response is zero due
to the compensating non-zero magnetoelectric response of each magnetic
sublattice. This calculation also highlights that the the orbital contribution
to the response is remarkably larger than the spin one and that each
sublattice has a rather large total magnetoelectric response of 85 ps/m.
Furthermore, we provide an intuitive recipe to visualize the dynamical
magnetic effective charge, allowing to examine its multipolar nature which we
confirm by means of ab initio calculations. Given the large value of the local
response, we investigate the ferromagnetic phase as well, which gives a giant
magnetoelectric response of about 1000 ps/m and coming mainly from the spin
contribution this time. Finally, we discuss the possible reasons for such a
large magnetoelectric response in BiCoO3 and propose possible strategies to
unveil this potentially large response.

##  I Introduction

Magnetoelectricity, that is the coupling between the electric and magnetic
orders in materials, experienced an increased interest in the past 20 years
[1, 2, 3] due to its promising technological applications [4, 5, 6, 7].
Although the idea of a magnetoelectric (ME) effect was first proposed by
Röntgen in 1888 [8], with Pierre Curie having emphasizing in 1894 by symmetry
analysis that intrinsic MEs could exist [9]; the name magnetoelectricity was
coined only a few decades later by Debye [10]. Cr2O3 was the first single
crystal to be theoretically proposed [11] and, soon after, experimentally
demonstrated [12, 13, 14, 15] to show a ME effect.

ME effects show up as cross-coupling terms in a Landau free energy expansion
in terms of applied magnetic ($\boldsymbol{\mathcal{H}}$) and electric
($\boldsymbol{\mathcal{E}}$) fields:

$F_{\text{ME}}=-\alpha_{ij}\mathcal{E}_{i}\mathcal{H}_{j}-\frac{1}{2}\beta_{ijk}\mathcal{E}_{i}\mathcal{E}_{j}\mathcal{H}_{k}-\frac{1}{2}\gamma_{ijk}\mathcal{H}_{i}\mathcal{H}_{j}\mathcal{E}_{k}+\dots.$
(1)

Particularly, here the first term, identified by the second-rank tensor
$\alpha_{ij}$ defined as

$\alpha_{ij}=\mu_{0}\frac{\delta
M_{j}}{\delta\mathcal{E}_{i}}\biggr{\rvert}_{\mathcal{H}}=\frac{\delta
P_{i}}{\delta\mathcal{H}_{j}}\biggr{\rvert}_{\mathcal{E}},$ (2)

describes the so-called linear ME effect, which results in a net magnetization
$\mathbf{M}$ (polarization $\mathbf{P}$) linear in the applied electric
(magnetic) field’s amplitude. On the other hand, $\beta_{ijk}$ and
$\gamma_{ijk}$ in Eq. (1) are third-rank tensors accounting for the (non-
linear) quadratic ME coupling. Because of the following upper bound for the ME
response [16]:

$\alpha_{ij}^{2}\leq\epsilon_{0}\mu_{0}\epsilon_{ii}\mu_{jj},$ (3)

where $\epsilon_{0}$, $\mu_{0}$, $\epsilon_{ii}$ and $\mu_{jj}$ are the vacuum
and material permittivities and permeabilities respectively, the ME coupling
is expected to be large in multiferroic materials due to the combined large
electric permittivity and large magnetic permeability [17, 18, 19, 20, 21,
22]. For this reason a great effort has been made to understand the different
mechanisms of magnetic and ME interactions in multiferroics from theory,
simulations and experiments [23, 24, 25, 26, 27, 28].

On the computational side, it is convenient to split the full ME response into
electronic, lattice-mediated and strain-mediated contributions
$\alpha^{tot}=\alpha^{\text{elec}}+\alpha^{\text{latt}}+\alpha^{\text{str}}$.
Furthermore, each term is usually separated into the contributions coming from
the spin and orbital degrees of freedom. Such calculations have shown that the
lattice mediated terms can have a significant contribution to the full
response when close to a ferroelectric phase transition where a polar soft
mode is at play [24, 29, 20] and that the electronic contribution is large
when close to a metal-insulator phase transition [20]. More importantly, when
the quenching of the orbital angular momentum is absent, a strong contribution
from the orbital magnetization can also be at play [30, 31, 23, 32]. For
example, this is the case for LiFePO4 where the lattice-mediated orbital
contribution to $\alpha_{\parallel}$ accounts almost entirely for the in-plane
ME response at T$=$0 K [30].

Over the past couple of decades, an increasing interest to a convenient
alternative description of linear ME responses based on microscopic markers
named ME multipoles has been observed in the literature [33, 34, 35, 36].
These markers identify the second-order contribution in a multipole expansion
of the energy of a magnetization density $\boldsymbol{\mu}(\mathbf{r})$
interacting with a non-uniform magnetic field
$\boldsymbol{\mathcal{H}}(\mathbf{r})$:

$\displaystyle E$ $\displaystyle=$
$\displaystyle-\int\boldsymbol{\mu}(\mathbf{r})\cdot\boldsymbol{\mathcal{H}}(\mathbf{r})\,d^{3}\mathbf{r}$
$\displaystyle=$
$\displaystyle-\boldsymbol{\mathcal{H}}(0)\cdot\int\boldsymbol{\mu}(\mathbf{r})\,d^{3}\mathbf{r}-\partial_{i}\mathcal{H}_{j}(0)\int
r_{i}\mu_{j}(\mathbf{r})\,d^{3}\mathbf{r}\quad-\ldots$

The tensor $\mathcal{M}_{ij}=\int r_{i}\mu_{j}(\mathbf{r})$ appearing in Eq.
(LABEL:eq:dipole), is called ME multipole tensor and has a one-to-one link to
the linear ME response $\alpha_{ij}$, because they both break separately
spatial- ($\mathcal{I}$) and time-inversion ($\mathcal{T}$) symmetries, while
preserving the combined symmetry $\mathcal{TI}$.

ME multipoles can be defined both for the whole unit cell or for single atoms,
depending on whether the integration in Eq. (LABEL:eq:dipole) is performed
over the unit cell or over a sphere centered at the atomic site, respectively
[34]. For the purpose of the present work, we will only refer to the latter,
i.e., the so-called atomic-site multipoles, which identify the contribution to
the ME response by each single atom. Practically, a ferroic arrangement of a
multipole appearing in given entries $\mathcal{M}_{ij}$ of the multipole
tensor results in corresponding non-zero entries $\alpha_{ij}$ of the linear
ME tensor. For example, ferroically ordered ME monopoles $a$ identify a net
diagonal and isotropic linear ME response.

Very recently, it has been shown that a vanishing net ME response does not
prevent a local, atomic ME response [36]. In fact, there can be cases where
the local ME responses of the oppositely aligned antiferromagnetic sublattices
are opposite to each other, thus giving a zero bulk net response [36]. Such
effect, named anti-ME, corresponds microscopically to antiferroically ordered
ME multipoles.

In this work, we explore the anti-ME effect in the multiferroic material
BiCoO3 by means of ab initio density functional theory (DFT). We specifically
consider the electronic and lattice contributions to the ME coupling in the
C-type antiferromagnetic (C-AFM) and ferromagnetic (FM) phases of BiCoO3. This
work is organized as follows: in Sec. II we summarize the theoretical approach
and list the computational parameters for each code used in our calculations;
in Sec. III we describe the structural and magnetic properties of BiCoO3; in
Sec. IV we focus on the Born effective charges and dynamical magnetic
effective charges and their contribution to the ME tensor; in Sec. V we
rationalize the absence of ME response in the C-AFM ground-state state of
BiCoO3 and the potentially large ME response of the FM phase; in Sec. VI we
give a conclusion and we provide some suggestions on how the large ME response
of the FM phase could be unlocked in BiCoO3.

##  II Methods

###  II.1 Theoretical formalism

####  II.1.1 Electronic and lattice-mediated ME responses

The electronic (i.e., clamped-ions) contribution to the ME tensor was
evaluated by applying a Zeeman magnetic field and computing the induced
polarization, without further relaxing the ions. In this approach, the
polarization $\mathbf{P}$ induced by a magnetic field
$\boldsymbol{\mathcal{H}}$ is

$P_{i}(\boldsymbol{\mathcal{H}})=\frac{1}{\Omega}\alpha_{ij}^{\text{elec}}\mathcal{H}_{j},$
(5)

where $\alpha_{ij}^{\text{elec}}$ is the electronic contribution to the ME
response, and $\Omega$ is the unit-cell volume.

The lattice-mediated contribution to $\alpha$ can be obtained in principle
within the same Zeeman field approach by further relaxing the atomic positions
in presence of the field. However, we note that in this way the orbital
contribution to $\alpha_{ij}^{\text{elec}}$ would not be accounted for because
the Zeeman field acts only on the spin degrees of freedom. Thus, we used a
method related to the one proposed in Ref. [24] where we estimated the
lattice-mediated part of the ME response as follows:

$\alpha^{\text{latt}}_{ij}=\frac{\mu_{0}}{\Omega}\sum_{n=1}^{N_{\text{IR}}}\frac{S_{n,ij}}{\omega_{n}^{2}},$
(6)

with $\omega_{n}$ the eigenvalues of the dynamical matrix (i.e., phonon
angular frequencies 111Note that the phonon angular frequencies $\omega$ are
related to the frequencies $\nu$, commonly expressed in THz or cm-1, by
$\omega=2\pi\nu$.) corresponding to the $N_{\text{IR}}$ infrared-active phonon
modes, whereas $S_{n,ij}$ is, what we will call the ME mode-oscillator
strength tensor, defined as follows:

$S_{n,ij}=\left(\sum_{\kappa,i^{\prime}}Z^{*\text{m}}_{\kappa,ii^{\prime}}{u}_{n,\kappa,i^{\prime}}\right)\times\left(\sum_{\kappa,j^{\prime}}Z^{*\text{e}}_{\kappa,j^{\prime}j}{u}_{n,\kappa,j^{\prime}}\right),$
(7)

where $\kappa$ identifies the atom in the unit cell, ${u}_{n,\kappa,j}$ the
normalized eigendisplacement, i.e., the generalized eigenvector of the
interatomic force constants matrix (see Supplementary Material for a more
detailed discussion), $Z_{\kappa}^{*\text{e}}$ and $Z_{\kappa}^{*\text{m}}$
the Born effective charge (BEC) and dynamical magnetic effective charge (DMC),
respectively, defined as:

$Z^{*\text{e}}_{\kappa,ij}=\Omega\frac{\partial
P_{i}}{\partial\tau_{\kappa,j}}=\frac{\partial^{2}E}{\partial\mathcal{E}_{i}\partial\tau_{\kappa,j}},$
(8)

i.e., the net polarization along direction $i$ induced by a displacement
$\tau$ of atom $\kappa$ along direction $j$, and

$Z^{*\text{m}}_{\kappa,ij}=\Omega\frac{\partial
M_{i}}{\partial\tau_{\kappa,j}}=\frac{\partial^{2}E}{\partial\mathcal{H}_{i}\partial\tau_{\kappa,j}},$
(9)

which corresponds to the net magnetization along direction $i$, induced by a
displacement $\tau$ of atom $\kappa$ along direction $j$ [29]. Hence, this ME
mode oscillator strength formula is similar to the so-called mode oscillator
strength of polar modes [38] but here for modes that are both infrared active
and magnetically active (i.e. the phonon mode carries a net magnetization).

We note that the scheme proposed by Iñiguéz in Ref. [24] to compute the
lattice-mediated ME response is based on the so-called mode effective charges
and mode dynamical magnetization, $\bar{Z}^{*\text{e}}_{n}$ and
$\bar{Z}^{*\text{m}}_{n}$ respectively, of the $n$-th phonon mode. These
identify the net polarization and net magnetization linearly induced by a
phonon mode, and are expressed as follows:

$\bar{Z}^{*\text{e}}_{n,i}=\sum_{\kappa,j}Z^{*\text{e}}_{\kappa,ij}\tilde{u}_{n,\kappa,j},$
(10)
$\bar{Z}^{*\text{m}}_{n,i}=\sum_{\kappa,j}Z^{*\text{m}}_{\kappa,ij}\tilde{u}_{n,\kappa,j}.$
(11)

It is worthwhile noting that the quantities identified by Eqs. (10) and (11)
are different from those inside the round brackets in Eq. (7), because the
eigendisplacements $u_{n,\kappa,j}$ and $\tilde{u}_{n,\kappa,j}$ are
normalized differently. We refer the reader to the Supplementary Material for
a more detailed discussion.

####  II.1.2 ME multipoles

The ME multipole tensor $\mathcal{M}$ introduced above is usually conveniently
decomposed into its spherical components [34], namely the ME monopole

$a=\frac{1}{3}{\cal
M}_{ii}=\frac{1}{3}\int\mathbf{r}\\!\cdot\boldsymbol{\mu}(\mathbf{r})\,d^{3}\mathbf{r},$
(12)

the toroidal moment vector

$\mathbf{t}=\frac{1}{2}\int\mathbf{r}\\!\times\boldsymbol{\mu}(\mathbf{r})\,d^{3}\mathbf{r},$
(13)

having components

$t_{i}=\frac{1}{2}\varepsilon_{ijk}{\cal M}_{jk}\quad,$ (14)

and the traceless quadrupole tensor

$\displaystyle q_{ij}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\left({\cal
M}_{ij}+{\cal M}_{ji}-\frac{2}{3}\delta_{ij}{\cal M}_{kk}\right)$ (15)
$\displaystyle=$
$\displaystyle\frac{1}{2}\int\left[r_{i}\mu_{j}(\mathbf{r})+r_{j}\mu_{i}(\mathbf{r})-\frac{2}{3}\delta_{ij}\mathbf{r}\\!\cdot\boldsymbol{\mu}(\mathbf{r})\right]d^{3}\mathbf{r},$

such that the overall ${\cal M}$ tensor is:

$\begin{pmatrix}a+\frac{1}{2}q_{x^{2}-y^{2}}-\frac{1}{2}q_{z^{2}}&t_{z}+q_{xy}&t_{y}+q_{xz}\\\
-t_{z}+q_{xy}&a-\frac{1}{2}q_{x^{2}-y^{2}}-\frac{1}{2}q_{z^{2}}&-t_{x}+q_{yz}\\\
-t_{y}+q_{xz}&t_{x}+q_{yz}&a+q_{z^{2}}\end{pmatrix}.$ (16)

As stated earlier, ferroically ordered ME multipoles appearing in a given
entry $ij$ of the multipole tensor $\mathcal{M}$ result in a net ME response
$\alpha_{ij}$. For example, ferroically ordered monopoles $a$ imply a net
diagonal and isotropic ME response.

###  II.2 Ab initio calculations

First-principles calculations were performed within DFT and density functional
perturbation theory (DFPT) [39], as implemented in the Vienna Ab Initio
Simulation package (VASP 5.4.4) [40, 41] and ABINIT (9.11.0) [42] codes. We
crosschecked between VASP and ABINIT codes to compute the lattice-mediated
contribution to the ME response (Eq. (6)). Further on, we used ABINIT to
compute the electronic contribution to the ME response (Eq. (5)). In both
codes, phonon eigenvalues and eigenvectors at zone center and BECs were
computed using DFPT, while DMCs were computed using a finite differences (FD)
approach. The calculation of orbital magnetization was carried out using the
muffin-tin approximation by summing the partial atomic contributions in the
PAW spheres and neglecting the contribution in the interstitial region. We
note that to get a complete evaluation of the orbital magnetization, including
the interstitial contributions, the modern theory of magnetization would be
required [43, 44, 45].

Next, we present the computational parameters used for both codes.

####  II.2.1 VASP

Lattice relaxation calculations were done within the generalized gradient
approximation (GGA) with a revised exchange–correlation functional for solids
(PBEsol) [46] and a relaxation threshold of $10^{-8}$ eV/Å on the ionic forces
was used. Electrons were described using Projector Augmented-Wave (PAW)
pseudopotentials, with the following valence electrons configurations: Bi
$6s^{2}6p^{3}$ (dataset Bi), Co $4s^{1}3d^{8}$ (dataset Co_sv), O
$2s^{2}2p^{4}$ (dataset O). The pseudo-wave functions were expanded into plane
waves [47] with a kinetic energy cut-off of 650 eV. The Brillouin zone (BZ) of
the magnetic unit cell, corresponding to a $\sqrt{2}$ ×$\sqrt{2}$×$1$ crystal
supercell (see Section III, was sampled by a 12$\times$12$\times$12 Monkhorst-
Pack k-point mesh [48].

DMCs were computed with the FD technique explained above, with self-consistent
field (SCF) calculations performed at the distorted structures obtained by
moving the atoms with a set of finite displacements of $\pm$ 0.01, 0.03, and
0.05 Å for both the in-plane and out-of-plane responses. The non-collinear SCF
calculations with spin-orbit coupling (SOC) to compute the DMCs were performed
within the local-density approximation (LDA), using into account the Slater
exchange and Perdew-Zunger parametrization of Ceperley-Alder Monte-Carlo
correlation data [49]. We used a strict convergence threshold of $10^{-12}$ eV
of the total energy on the self-consistent cycle and added 500-2000 steps to
allow for the convergence of the magnetic moment in order to determine the
DMCs accurately. Correlation effects were dealt with by applying a Hubbard
correction [50] to the d orbitals of Cobalt, with a U parameter of 6 eV, which
opens a band gap (1.42 eV) in agreement with previously reported theoretical
values [51, 52] and comparable to experimental values from X-ray emission and
absorption spectroscopy (XES/XAS) [53] (1.7 eV).

The ME and higher-order multipoles were computed by decomposing the DFT+U
occupation matrix into spherical tensor moments [54], as obtained from an
extension of the Multipyles Python package 222See
https://github.com/materialstheory/multipyles..

##### Visualization of DMC-related perturbed magnetization densities.

In order to create an isosurface to illustrate the DMCs, we calculated at each
point of the $96\times 96\times 80$ Fast Fourier Transformation (FFT) grid the
difference in magnetization density between the perturbed and ground states.
The display limit value, called the isosurface level, was set in order to get
a similar surface area for both the linear and quadratic contributions shown
in Fig. 2, discussed below in Sec. IV.2, for visual clarity. In order to
decompose the total response into linear and quadratic contributions (Fig. 2),
we first worked out the average of the magnetization density from two opposite
atomic displacements, which allowed us to get the quadratic part. Then, we
subtracted it from the total induced magnetization density to get the linear
contribution.

####  II.2.2 ABINIT

Calculations with ABINIT were performed in the LSDA+U framework [56] using the
JTH PAW atomic potentials (v1.1) [57], with SOC included and starting from the
relaxed structures obtained in VASP. A Hubbard-U correction of 5.0 eV was
applied to the 3d Co shell, while although different from the one used in
VASP, it gives the same band gap (since the effective U is affected by the
radii of the PAW spheres, which are different in the two codes). The Kohn-Sham
pseudo-wave functions were expanded over a plane-wave basis set with a kinetic
energy cut-off of 20 Ha (544 eV). We used a strict convergence threshold of
$2.7\times 10^{-12}$ eV for the total energy on the self-consistent cycle, to
determine the DMCs accurately. The BZ of the magnetic unit cell was sampled
with a Monkhorst-Pack k-point grid of 8x8x8 points. For DFPT calculations
[38], we fully re-optimized the structure, by using a tolerance of 1.3$\times
10^{-6}$ GPa on the internal pressure and of 3$\times 10^{-7}$ eV/Å on the
forces.

##  III Structural and magnetic properties

BiCoO3 crystallizes in a supertetragonal phase described by the space group
$P4mm$ (point group 4mm) with a much more pronounced tetragonality ($c/a$ =
1.267) than the reference PbTiO3 crystal for this phase. The crystal unit cell
contains 5 atoms, the Co atoms are situated in a tetragonal prism having four
O atoms at the corner of the base and one at the apex, side-sharing the base
with an opposing trigonal prism with the Bi atom at the apex. Particularly, Bi
atoms sit at the 1a (0, 0, 0) Wyckoff site, Co atoms are located at the 1b
Wyckoff position ($\frac{1}{2}$, $\frac{1}{2}$, zCo), whereas the apical and
basal oxygens are found at the 1b ($\frac{1}{2}$, $\frac{1}{2}$,
z${}_{O_{api}}$) and 2c ($\frac{1}{2}$, 0, z${}_{O_{bas}}$), respectively.

Figure 1: Schematic representation of a) the 5-atoms primitive cell of BiCoO3
without magnetism and b) the $\sqrt{2}$×$\sqrt{2}$×1 C-AFM unit cell. Green
arrows show the direction of the Co magnetic moment.

The pronounced tetragonality, typical of the super-tetragonal structure,
allows for a giant polarization of 162.1 $\mu$C/cm2 as estimated from the BECs
calculated with VASP and using atomic displacements relative to a reference
structure with space group $P4/mmm$, i.e., the closest centrosymmetric
structure. This calculated value is in good agreement with the previously
reported value of 179 $\mu$C/cm² [52], obtained with ab initio Berry-phase
calculations.

BiCoO3 also undergoes an antiferromagnetic (AFM) phase transition at the Néel
temperature $T_{N}=470~{}K$ [58]. The magnetic ground state of BiCoO3 has a
propagation vector $k=(1/2,1/2,0)$, which leads to a C-AFM phase described by
a magnetic unit cell with 10 atoms, obtained as a
$\sqrt{2}\times\sqrt{2}\times 1$ replica of the crystal unit cell described
earlier (see Fig. 1). We performed structural relaxations, followed by self-
consistent total energy calculations at the relaxed structure for different
magnetic orders, namely A-AFM, C-AFM, G-AFM, with propagation vector
$k=(0,0,1/2)$, $k=(1/2,1/2,0)$, $k=(1/2,1/2,1/2)$, and ferromagnetic (FM). Our
results, summarized in Tab. 1, show that the C-AFM is the most energetically
favorable magnetic state in agreement with previous experimental measurements
[58]. Furthermore, our calculations with SOC give that BiCoO3 has an easy axis
along the $z$ axis (the tetragonal axis), which we will assume through the
rest of the manuscript. The magnetic moment of Co atoms is found to be
3.16/3.24 $\mu_{\text{B}}$ (3.0/3.1 $\mu_{\text{B}}$ from spin and 0.16/0.14
$\mu_{\text{B}}$ from orbital moment) respectively with VASP and ABINIT.

Formally the d6 configuration leads to a magnetic moment of 4$\mu_{\text{B}}$,
but the strong hybridization between the Co d and O p states makes the
magnetic moment of Cobalt dilute on the apical site. This result is also seen
with muffin-tin spheres calculations [52]

Table 1: Energies and lattice parameters of the magnetic phases of BiCoO3 calculated with VASP and ABINIT (in italic font), compared to the experimental values (last row). All the reported values and space groups are for magnetic moments oriented along the $z$ direction. | Mag Space
---
group
| Magnetic
---
Phase
E(meV/f.u.) |  $a$ (Å) |  $c$ (Å) | $c/a$
$P_{p}4^{\prime}mm^{\prime}$ | C-AFM | 0 | 7.390 | 9.362 | 1.267
$C_{I}m^{\prime}m^{\prime}2$ | G-AFM | 17 | 7.389 | 9.360 | 1.267
$P_{2c}m^{\prime}m^{\prime}2$ | A-AFM | 162 | 7.432 | 9.408 | 1.266
$P4m^{\prime}m^{\prime}$ | FM | 192 | 7.421 | 9.471 | 1.276
|  |  | 7.350 | 9.400 | 1.279
Exp [58] | C-AFM |  | 7.459 | 9.448 | 1.267

In Tab. 1 we also report the magnetic space groups (in the OG setting) for all
magnetic phases, assuming the Co magnetic moment aligned along the $z$
direction. Among the four magnetic states analyzed here, only the FM state
allows a linear ME response, which is diagonal by symmetry, i.e.
$\alpha_{xx}=\alpha_{yy}\neq\alpha_{zz}$. On the other hand, time-reversal
combined with a fractional translation implies that a net linear ME response
is forbidden by symmetry in any of the AFM states. To understand this more
precisely at the microscopic level, next we analyze the ME response (or the
lack thereof) of BiCoO3 in both the FM and the C-AFM phases.

##  IV Electronic and lattice responses

As we have discussed in Sec. II.1, to probe the ME response, we can combine
DFPT and finite difference techniques (i.e. by determining the BEC and DMC,
and using Eq. (6)) or an applied magnetic field technique and extract the
induced polarization versus the field amplitude. In the following, we apply
both of them to BiCoO3 in its C-AFM and FM magnetic phases and we compare the
results. We start from the calculations of the BECs and DMCs.

###  IV.1 Born effective charges

Table 2: Calculated Born effective charge tensor components of BiCoO3, obtained with both VASP and ABINIT (in italic font).  Atom | Wyc. | Born Effective Charges (e)
---|---|---
|  |  | $\partial/\partial\mathcal{E}_{x}$ | $\partial/\partial\mathcal{E}_{y}$ | $\partial/\partial\mathcal{E}_{z}$
Bi | 1a | $\partial/\partial\tau_{x}$ | 5.08 | (5.07) | 0 |  | 0 |
$\partial/\partial\tau_{y}$ | 0 |  | 5.08 | (5.07) | 0 |
$\partial/\partial\tau_{x}$ | 0 |  | 0 |  | 3.81 | (3.72)
Co | 1b | $\partial/\partial\tau_{x}$ | 2.63 | (2.30) | 0 |  | 0 |
$\partial/\partial\tau_{y}$ | 0 |  | 2.63 | (2.30) | 0 |
$\partial/\partial\tau_{x}$ | 0 |  | 0 |  | 3.60 | (3.43)
O${}_{\text{api}}$ | 1b | $\partial/\partial\tau_{x}$ | -2.22 | (-2.13) | 0 |  | 0 |
$\partial/\partial\tau_{y}$ | 0 |  | -2.22 | (-2.13) | 0 |
$\partial/\partial\tau_{x}$ | 0 |  | 0 |  | -3.26 | (-3.19)
O${}_{\text{bas}}$ | 2c | $\partial/\partial\tau_{x}$ | -2.75 | (-2.62) | 0.02 | (-0.14) | 0 |
$\partial/\partial\tau_{y}$ | 0.02 | (-0.14) | -2.75 | (-2.62) | 0 |
$\partial/\partial\tau_{x}$ | 0 |  | 0 |  | -2.07 | (-1.98)

Our calculated BECs are reported in Tab. 2 for both VASP and ABINIT codes.
There is a substantial agreement between the two codes with maximum
discrepancies of 15% of differences on few components, which might be due to
the different radii of the PAW spheres or other atomic data parameters
difference. We can also observe that all the BEC are diagonal (only the
O${}_{\text{bas}}$ has small off-diagonal terms along the $x$ and $y$
directions), i.e., the charge rearrangement due to atomic displacements takes
place along the direction of the atomic displacement. Further on, BECs allow
us to compute the mode effective charge Eq. (10) of polar phonons, which enter
in Eq. (6).

In the non-centrosymmetric space group $P4mm$, the infrared (IR) and Raman
active modes coincide, i.e., a mode is silent or both Raman and IR active. The
12 zone-center optical phonon modes of the 5-atoms unit cell are classified by
symmetry according to the irreducible representations of $P4mm$ in the
following way: 3$\Gamma_{1}$+$\Gamma_{3}$+8$\Gamma_{5}$. $\Gamma_{1}$
($A_{1}$) modes are IR active along the $z$ direction and non-degenerate,
whereas $\Gamma_{5}$ ($E$) modes are IR active in the $xy$ plane and two-fold
degenerate. The $\Gamma_{3}$ ($B_{1}$) mode is silent. In Tab. 3 we report the
frequency and the calculated mode effective charge for the IR active modes.
The values obtained with VASP and ABINIT are consistent with one another.

Table 3: Frequency and mode effective charges of IR active modes of BiCoO3, as calculated from VASP (regular font) and ABINIT (italic font). $\nu$(cm-1) | Character | Mode effective charge (e)
---|---|---
$x$ | $y$ | $z$
83 | 65 | $\Gamma_{5}$ | 4.91 | 4.79 | 0 | 0 | 0 | 0
83 | 65 | $\Gamma_{5}$ | 0 | 0 | 4.91 | 4.79 | 0 | 0
218 | 218 | $\Gamma_{1}$ | 0 | 0 | 0 | 0 | -2.79 | -2.71
220 | 216 | $\Gamma_{5}$ | -2.10 | -1.50 | 0 | 0 | 0 | 0
220 | 216 | $\Gamma_{5}$ | 0 | 0 | -2.10 | -1.50 | 0 | 0
402 | 404 | $\Gamma_{5}$ | -5.61 | -5.41 | 0 | 0 | 0 | 0
402 | 404 | $\Gamma_{5}$ | 0 | 0 | -5.61 | -5.41 | 0 | 0
417 | 430 | $\Gamma_{1}$ | 0 | 0 | 0 | 0 | 5.84 | 5.37
546 | 542 | $\Gamma_{5}$ | -1.95 | -1.70 | 0 | 0 | 0 | 0
546 | 542 | $\Gamma_{5}$ | 0 | 0 | -1.95 | -1.70 | 0 | 0
672 | 699 | $\Gamma_{1}$ | 0 | 0 | 0 | 0 | -4.23 | -4.39
Figure 2:  (a) Plot of the induced spin magnetization along the $z$ direction versus the amplitude of the Co atom displacement along the $z$ direction. Only the induced magnetization along the $z$ direction is non-zero and it shows a non-linear behavior on top of the linear response (highlighted by the dashed black line. Green and red arrows identify the linear and the non-linear contributions to the induced magnetization, respectively. Panels b, c, and d, give a visualization of the total, linear and non-linear induced magnetization density, respectively. Table 4: Linear spin ($\vec{S}$) and orbital ($\vec{L}$) dynamical magnetic effective charges (10${}^{-2}\mu_{\text{B}}/\text{\r{A}}$) for each inequivalent Wyckoff sites of BiCoO3 as calculated with VASP and ABINIT (ABINIT values are in brackets and in italic). $\mathcal{H}_{i}$ refers to the magnetic field in direction $i$ and $\tau_{j}$ to the atom displacement along the direction $j$. Atom | Wyckoff | $\vec{S}$ |  $\vec{L}$
---|---|---|---
|  |  | $\partial/\partial\mathcal{H}_{x}$ | $\partial/\partial\mathcal{H}_{y}$ | $\partial/\partial\mathcal{H}_{z}$ |  | $\partial/\partial\mathcal{H}_{x}$ | $\partial/\partial\mathcal{H}_{y}$ | $\partial/\partial\mathcal{H}_{z}$
Bi | 1a | $\partial/\partial\tau_{x}$ | 6.93 | (6.93) | 0 |  | 0 |  | $\partial/\partial\tau_{x}$ | 9.83 | (8.32) | 0 |  | 0 |
$\partial/\partial\tau_{y}$ | 0 |  | -6.93 | (-6.93) | 0 |  | $\partial/\partial\tau_{y}$ | 0 |  | -9.83 | (-8.32) | 0 |
$\partial/\partial\tau_{z}$ | 0 |  | 0 |  | 0 |  | $\partial/\partial\tau_{z}$ | 0 |  | 0 |  | 0 |
Co | 1b | $\partial/\partial\tau_{x}$ | 7.25 | (8.49) | 0 |  | 0 |  | $\partial/\partial\tau_{x}$ | 24.63 | (22.57) | 0 |  | 0 |
$\partial/\partial\tau_{y}$ | 0 |  | 7.25 | (8.49) | 0 |  | $\partial/\partial\tau_{y}$ | 0 |  | 24.63 | (22.57) | 0 |
$\partial/\partial\tau_{z}$ | 0 |  | 0 |  | -0.91 | (-0.67) | $\partial/\partial\tau_{z}$ | 0 |  | 0 |  | -15.01 | (-15.34)
O${}_{\text{api}}$ | 1b | $\partial/\partial\tau_{x}$ | 2.80 | (1.74) | 0 |  | 0 |  | $\partial/\partial\tau_{x}$ | -6.12 | (-5.52) | 0 |  | 0 |
$\partial/\partial\tau_{y}$ | 0 |  | 2.80 | (1.74) | 0 |  | $\partial/\partial\tau_{y}$ | 0 |  | -6.12 | (-5.52) | 0 |
$\partial/\partial\tau_{z}$ | 0 |  | 0 |  | 0.11 | (0.10) | $\partial/\partial\tau_{z}$ | 0 |  | 0 |  | -0.57 | (-1.83)
O${}_{\text{bas}}$ | 2c | $\partial/\partial\tau_{x}$ | 0 |  | 0 |  | 1.1 | (0.4) | $\partial/\partial\tau_{x}$ | 0 |  | 0 |  | 13 | (6.1)
$\partial/\partial\tau_{y}$ | 0 |  | 0 |  | -1.1 | (-0.4) | $\partial/\partial\tau_{y}$ | 0 |  | 0 |  | -13 | (-6.1)
$\partial/\partial\tau_{z}$ | -5.8 | (-5.9) | 5.8 | (5.9) | 0 |  | $\partial/\partial\tau_{z}$ | -3.3 | (-3.8) | 3.3 | (3.8) | 0 |

###  IV.2 Dynamical magnetic effective charges

We calculated the DMC from finite differences by displacing the atoms along
each symmetry-inequivalent directions (See Supplementary Material Fig. S1. for
the raw data). The range of amplitudes chosen allows us to capture linear as
well as higher-order, non-linear, contributions to the induced magnetization
as exemplified in Fig. 2 (a) for the case of Co atom.

We establish that most of the responses are linear, but a few of them are
significantly non-linear. This is the case for the change in orbital
magnetization along the $z$ direction ($\Delta m^{\text{orb}}_{z}$) for Co
moving in-plane and the in-plane change of orbital magnetization ($\Delta
m^{\text{orb}}_{x}$ and $\Delta m^{\text{orb}}_{y}$) for O${}_{\text{bas}}$
moving in the in-plane direction (see Supplementary Material Fig. S1).

Figure 3: Schematic representation of the DMC pseudotensor as viewed from
[001] (left panel) and [100] (right panel) directions. The color of the
vectors surrounding Bi, Co and apical Oxygen atom represent the in-plane
direction of the induced magnetization when the atom is displaced in a given
direction that cross a sphere of arbitrary radius, such that we end up with a
3D vector field for each possible displacement direction around each atom. The
vectors surrounding the Basal oxygen are blue when the out-of-plane component
of the induced magnetization is along z and yellow for the opposite direction.

In Table 4 we summarize the spin and orbital contributions to the DMC tensor
components, as obtained in both VASP and ABINIT codes from the linear term of
the fitted induced magnetization as a function of the atomic displacement. Our
results show that the orbital contribution is larger than the spin
contribution for most of the DMC components. This effect is not common and it
is particularly enhanced for Co, which can be explained because the Co3+
cation with d6 orbital filling is known to potentially have an unquenched
orbital moment in crystals. Furthermore, we note that the computed DMCs,
especially for the orbital contribution, are remarkably larger than those of
the paradigmatic ME material, Cr2O3. In fact, the largest orbital DMC in
BiCoO3 is approximately $0.25\,\mu_{\text{B}}/\text{\r{A}}$ for Co, two orders
of magnitude larger than the largest orbital DMC reported for Cr2O3
($0.0025\,\mu_{\text{B}}/\text{\r{A}}$ for Cr) [29].

As explained in Ref. [35], the DMCs are proportional to the local linear ME
response. For this reason, the shape of the DMC for each atom can be
conveniently explained by either (i) the correspondence to the atomic-site ME
multipoles introduced earlier in Sec. I, or similarly (ii) the shape of the
local linear ME response, i.e., the linear ME tensor relative to the point
symmetry of the Wyckoff site the atom occupies. We clarify this point by
taking, as an example, the case of the atoms of C-AFM BiCoO3. In Table 5 we
report the computed ME multipoles for BiCoO3. Bi atoms display only a non-zero
$x^{2}-y^{2}$ ME quadrupole, which results in a diagonal DMC with only one
independent parameter, namely $Z^{*\text{m}}_{yy}=-Z^{*\text{m}}_{xx}$
(whereas $Z^{*\text{m}}_{zz}=0$). This corresponds to the shape of the linear
ME tensor for Bi’s Wyckoff site 2$a$, with local symmetry
$4^{\prime}mm^{\prime}$. Co and the apical O both have the same allowed ME
multipoles, namely the ME monopole $a$ and the $z^{2}$ quadrupole, because
they sit in the same Wyckoff site 2$b$. Their allowed multipoles, according to
Eq. (16), explain the diagonal shape of the DMC, with two independent
parameters $Z^{*\text{m}}_{yy}=Z^{*\text{m}}_{xx}$ and $Z^{*\text{m}}_{zz}\neq
0$, in agreement with the linear ME tensor for the point group
$4m^{\prime}m^{\prime}$ (i.e., the site symmetry of Wyckoff position 2$b$).
Finally, basal oxygens have non-vanishing toroidal moments $t_{x}$ and
$t_{y}$, as well as $xz$ and $yz$ quadrupoles, which result in the DMC with
the shape reported in Table 4. This has more independent components than the
other atoms because the site symmetry of basal O (Wyckoff site 4$c$) is lower
($2^{\prime}.m^{\prime}m$).

We remark that DMCs are anisotropic tensors, implying that their shape depends
on the orientation of the reference framework. Practically, if the axes are
transformed by a rotation $\mathcal{R}$, the DMCs transform accordingly as
$Z^{\prime*\text{m}}=\mathcal{R}^{-1}Z^{*\text{m}}\mathcal{R}$.

Table 5: Calculated atomic-site ME multipoles (in units of $10^{-2}\,\mu_{\text{B}}\text{\AA}$) for BiCoO3. Atom | Monopole | Toroidal | Quadrupole
---|---|---|---
|  | moments | moments
| $a$ | $t_{x}$ | $t_{y}$ | $t_{z}$ | $\mathcal{Q}_{x^{2}-y^{2}}$ | $\mathcal{Q}_{z^{2}}$ | $\mathcal{Q}_{xy}$ | $\mathcal{Q}_{xz}$ | $\mathcal{Q}_{yz}$
Bi | – | – | – | – | $5.8\times 10^{-3}$ | – | – | – | –
Co | 1.1 | – | – | – | – | 1.0 | – | – | –
O${}_{\text{api}}$ | \- 1.6 | – | – | – | – | -1.7 | – | – | –
O${}_{\text{bas}}$ | – | 1.8 | 1.8 | – | – | – | – | 2.3 | -2.3

In Fig. 3 we propose a schematic way to visualize the DMCs. Particularly, for
a given set of points r centered around each atom, with r defined relative to
the equilibrium position of the atom, we show with an arrow the magnetization
induced by displacing the atom by r and its color is function of its
direction. Clearly, according to the definition of DMC given in Eq. (9), the
induced magnetization is $\Delta M_{i}=Z^{*\text{m}}_{ij}r_{j}$, proportional
to the DMC, thus the vector fields shown in Fig. 3 provide a visual
representation of the DMCs. From Fig. 3, we easily see that Co atom shows a
radial and isotropic (similar to an s atomic orbital) magnetization in the
$xy$ (i.e., $ab$ in the direct lattice coordinates) plane when shifted from
its equilibrium position. Furthermore, it shows a magnetization pattern
flattened along the $z$ (i.e., c) axis. Overall, the vector field representing
the induced magnetization density has the symmetry of a combined $s$ and
$d_{z^{2}}$ orbital, in agreement with Co showing a spherically symmetric
monopole $a$ and a $z^{2}$ quadrupole, as reported in Table 5. Apical O atoms
show a similar behavior to Co atom, to which they form robust bonds
characterized by strong hybridization. This hybridization causes a magnetic
moment transfer from Co atom to O atoms, contributing to the observed
magnetization pattern. In contrast to Co and apical O, Bi’s displacement
induces a magnetization with a shape resembling that of a $d_{x^{2}-y^{2}}$
atomic orbital, exhibiting lobes with opposite magnetization - one positive
and the other negative. This is in agreement with the reported $x^{2}-y^{2}$
ME quadrupole on Bi in Table 5. Finally, basal O induces a chiral/spiral
magnetization pattern, justified by its ME toroidal moment.

It is worthwhile noting that the results discussed so far are relative to one
magnetic sublattice (the C-AFM contains two magnetic sublattices of opposite
direction). The DMCs (ME multipoles) for the opposite magnetic sublattice are
equal in magnitude but opposite in sign to those reported in Table 4 (Table
5). This happens because the two sublattices are mapped into one another by
time-reversal symmetry ($\mathcal{T}$ combined with the $(1/2,1/2,0)$
fractional translation is a symmetry of the magnetic lattice), and the DMCs
and the ME multipoles change sign under such symmetry. On the other hand, the
BECs discussed in the previous subsection do not change sign in the two
sublattices because they are even under $\mathcal{T}$.

Having calculated the DMC for each atom, next we use Eq. (11) to compute the
mode DMCs. We focus particularly on the phonon modes at zone center, because
they are those contributing to the ME response. First, we note that the
10-atoms C-AFM magnetic unit cell has 30 phonon modes at $\Gamma$, which can
be grouped in two sets of 15 modes each. The first group contains modes found
at the zone center in the 5-atom (i.e., non magnetic or ferromagnetic) unit
cell, and thus characterized by in-phase atom displacements for the two
magnetic sublattices. Among the members of this group there are also IR
active, i.e., polar, phonon modes, which show a non-zero mode effective charge
as we discussed earlier in Table 3. The second set, instead, contains modes at
the zone-boundary high-symmetry point $M$ ($\textbf{k}=(1/2,1/2,0)$) of the
5-atom unit cell, which are folded to zone-center in the BZ of the magnetic
$\sqrt{2}\times\sqrt{2}\times 1$ supercell. Being zone-boundary modes of the
non-magnetic cell, they are characterized by atomic displacements with
opposing phases in the two magnetic sublattice cells. As such, these modes are
non-polar and their mode effective charges is null. As discussed earlier, each
atom in a given magnetic sublattice has an opposite DMC to the corresponding
atom in the other sublattice. Thus, these opposite contributions cancel each
other out in phonon modes of the first group, while they add up in phonon
modes of the second group, due to the opposite displacements of atoms
belonging to the two sublattices, which cancel the inversion of DMC. In Table
6 we report the mode dynamical magnetic effective charges for the magnetically
active modes, i.e., those at the M zone-boundary point of the non-magnetic
cell. In Table 6 we can see that the $M_{1}$ mode has a nonzero mode magnetic
effective charge along the $z$ direction and that the $M_{5}$ modes give
nonzero values in the $xy$ plane. All the other modes ($M_{2}$, $M_{3}$ and
$M_{4}$) are nill. Note that, because they have a degeneracy of two, the
projection of the pair of degenerate $M_{5}$ modes depends on the choice of
the eigenvector direction of the degenerate subspace, so the amplitude of the
dynamical magnetic effective charge must be considered as the sum of the pair.

Overall, the discussion above implies that polar phonon modes are not
magnetically active. This is true, in general, for any AFM with magnetic
symmetry described by a type-IV magnetic space group.

Table 6: Spin, orbital, and total mode DMCs for phonon modes at the M point of the BZ of the non-magnetic unit cell, which fold to $\Gamma$ in the BZ of the magnetic cell. Values obtained with VASP and ABINIT (italic font). $\nu$ (cm-1) | Character | Spin | Orbital | Total
---|---|---|---|---
92 | 82 | M5 (x2) | 6.8 | 0.4 | 8.6 | 5.3 | 15.4 | 5.7
146 | 129 | M5 (x2) | 6.0 | 3.8 | 13.4 | 5.4 | 19.4 | 9.2
221 | 224 | M1 | 0.4 | 1.1 | 11.0 | 21.9 | 11.4 | 23.0
309 | 308 | M5 (x2) | 1.2 | 11.7 | 28.1 | 31.1 | 32.3 | 42.8
347 | 349 | M5 (x2) | 2.8 | 7.8 | 7.3 | 11.3 | 10.1 | 19.1
633 | 623 | M1 | 2.9 | 1.5 | 34.6 | 23.4 | 37.5 | 24.9
654 | 681 | M1 | 1.5 | 0.1 | 16.2 | 0.1 | 17.7 | 0.2
Table 7: Allowed magnetic octupoles and corresponding allowed entries of the D2MC tensor for BiCoO3. Note: in $Z^{*m(2)}_{ijk}$ symmetry under exchange of $j$ and $k$ is implied, i.e., if $Z^{*m(2)}_{ijk}$ is allowed, $Z^{*m(2)}_{ikj}$ is allowed as well. Atom | Allowed octupoles | Allowed D2MC entries
---|---|---
Bi |  $\mathcal{O}_{2}$, $\mathcal{Q}^{\tau}_{xy}$ |  $Z^{*m(2)}_{xxz}=-Z^{*m(2)}_{yyz}$, $Z^{*m(2)}_{zxx}=-Z^{*m(2)}_{zyy}$
Co |  $\mathcal{O}_{0}$, $t^{\tau}_{z}$ |  $Z^{*m(2)}_{xxz}=Z^{*m(2)}_{yyz}$, $Z^{*m(2)}_{zxx}=Z^{*m(2)}_{zyy}$, $Z^{*m(2)}_{zzz}$
O${}_{\text{api}}$ |  $\mathcal{O}_{0}$, $t^{\tau}_{z}$ |  $Z^{*m(2)}_{xxz}=Z^{*m(2)}_{yyz}$, $Z^{*m(2)}_{zxx}=Z^{*m(2)}_{zyy}$, $Z^{*m(2)}_{zzz}$
O${}_{\text{bas}}$ |  $\mathcal{O}_{-3}=\mathcal{O}_{3}$, $\mathcal{O}_{-1}=-\mathcal{O}_{1}$, |  $Z^{*m(2)}_{xxx}=-Z^{*m(2)}_{yyy}$, $Z^{*m(2)}_{xxy}=-Z^{*m(2)}_{yxy}$, $Z^{*m(2)}_{yxx}=-Z^{*m(2)}_{xyy}$,
$\mathcal{Q}^{\tau}_{yz}=\mathcal{Q}^{\tau}_{xz}$ |  $Z^{*m(2)}_{xzz}=Z^{*m(2)}_{yzz}$, $Z^{*m(2)}_{zxz}=-Z^{*m(2)}_{zyz}$

###  IV.3 Non-linear Dynamical Magnetic Effective Charges

Fig. 2 unambiguously illustrates the feasibility of eliciting non-linear
responses where, for wide enough ranges of displacement amplitudes, the FD
approach is able to reliably capture higher-order contributions to the induced
magnetization. This is particularly true for cases where a linear induced
magnetization is forbidden and the quadratic one is the lowest-order response.
In most cases, the second-order contribution is much smaller than the linear
one; those where the second-order response is non-negligible and necessary to
precisely describe the induced magnetization are highlighted with dashed lines
in Fig. S1 of the Supplementary Material.

The quadratic induced magnetization can be quantified by the so-called second-
order dynamical magnetic effective charge (D2MC henceforth). This is a second-
order analogue of the DMC given in Eq. (9). It has been introduced in Ref.
[35] and it is defined as:

$Z^{*\text{m}(2)}_{ijk}=\Omega\frac{\partial^{2}M_{i}}{\partial\tau_{j}\partial\tau_{k}}.$
(17)

By construction, the D2MC is a third-rank pseudotensor, symmetric upon
exchange of $j$ and $k$. Given the definitions of the DMC and the D2MC, the
induced magnetization $\boldsymbol{\Delta M}$ as a function of the
displacement $\boldsymbol{\tau}$ reads

$\Delta
M_{i}=Z^{*\text{m}}_{ij}\tau_{j}+Z^{*\text{m}(2)}_{ijk}\tau_{j}\tau_{k},$ (18)

thus $Z^{*\text{m}(2)}$ can be obtained from a quadratic fit of
$\boldsymbol{\Delta M}(\boldsymbol{\tau})$. We note here that while
displacements along pure Cartesian directions allow to obtain the full DMC,
the same does not apply for the D2MC, and in fact only the entries with $j=k$
can be obtained in this way. To obtain the entries with $j\neq k$, linear
combinations of displacements along different Cartesian directions are
required.

While a quantitative evaluation of the full D2MC tensors is beyond the scope
of the present work, we remark that their allowed entries are determined by
symmetry and, similarly to the DMC discussed earlier, can be identified by a
link with specific multipoles of the magnetization density. As explained in
Ref. [35], the D2MCs are linked to the so-called magnetic octupoles, which
represent the next-order term in the expansion of Eq. (LABEL:eq:dipole) and
are identified by the third-rank pseudotensor $\mathcal{M}_{ijk}=\int
d^{3}\mathbf{r}\mu_{i}(\mathbf{r})r_{j}r_{k}$. Our computed octupoles in
BiCoO3 are reported in Table 7, following the naming convention of Ref. [35].
In the same Table, we indicate also the independent entries of the D2MCs
implied by the allowed octupoles.

As a final note, we remark that similarly to how the DMCs contribute to the
linear ME response (Eq. (6)), the D2MCs contribute to identify the second-
order ME response $\beta_{ijk}$, induced by an electric field, introduced in
Eq. (1). For a detailed discussion and the derivation we refer the interested
reader to the Supplementary Material of Ref. [36]. Although a quantitative
analysis of $\beta_{ijk}$ for each atom of BiCoO3 is not a target of this
work, it is important to note that, similarly to the linear ME response, the
second-order response is allowed locally, i.e., at the atomic level, but
vanishes at the unit cell level, due to the time-reversal symmetry combined
with a fractional translation connecting the two sublattices. In fact,
magnetic octupoles and, similarly, the D2MCs as well as $\beta_{ijk}$, change
sign under time-reversal, thus resulting in a second-order anti-ME response,
i.e., with equal magnitude but opposite sign in the two sublattices.

##  V Magnetoelectric response

Having calculated the BECs and the DMCs, we can compute the lattice
contribution to the ME response using Eq. (6). According to Eq. (6), to have a
large lattice-mediated ME response, a material needs to exhibit low-lying IR
phonon modes with correspondingly large mode effective charges and large
dynamical magnetic effective charges.

###  V.1 C-AFM phase

BiCoO3 cast large mode effective electric and magnetic charges but not on
matching phonon modes. Indeed, as reported above, the mode DMC is zero for the
polar modes, hence the total ME response is strictly zero. Let us look more
closely at why the ME response vanishes.

Figure 4: Change of polarization (along the three directions, i.e. $\Delta
P_{x}$, $\Delta P_{y}$ and $\Delta P_{z}$) versus applied magnetic field
amplitude along the $x$ direction as calculated in the C-AFM phase of BiCoO3
with ABINIT. Only the $\Delta P_{z}$ gives a response that is non-linear only.

At first sight it can be surprising that BiCoO3 hosts large atomic DMCs, which
result nonetheless in a null total ME response. This can be understood by
decomposing the total response into the response of each magnetic sublattice.
The time-reversal symmetry combined with the fractional translation
$(\frac{1}{2},\frac{1}{2},0)$ implies that the C-AFM order is described by a
propagation vector $\mathbf{q}=(\frac{1}{2},\frac{1}{2},0)$. The resulting
phase factor describing the magnetization density and, likewise, the induced
magnetization, implies that the DMCs are opposite for the two sublattices,
while the BECs have the same sign and magnitude. Thus, according to Eq. (6),
the two sublattices have non-vanishing, and opposite in sign ME responses,
resulting overall in a net vanishing ME response. This qualifies BiCoO3 as an
instance of the recently introduced anti-ME case [36]. We find that the
amplitude of the (lattice-mediated) ME response of each sublattice is
$\alpha_{xx}=85$ ps/m and $\alpha_{zz}=1$ ps/m. The $\alpha_{zz}$ response is
nonzero through the orbital contribution only, the spin response being zero as
expected [59, 30]. We remark that if the magnetic unit cell was the same as
the crystallographic one, the propagation vector would be $\mathbf{q}=(0,0,0)$
and hence the ME response would be nonzero.

In Fig. 4 we report the change of polarization of the C-AFM phase of BiCoO3
under an applied Zeeman magnetic field along the $x$ direction as obtained
with the ABINIT code. We confirm that the linear ME response is strictly zero
as no response is observed along the $x$ and $y$ directions ($\Delta P_{x}$
and $\Delta P_{y}$ in Fig. 4) and a pure quadratic response is observed along
the $z$ direction ($\Delta P_{z}$, blue curve in Fig. 4). Following Eq. (1),
the quadratic response observed along the $z$ direction corresponds to the
$\gamma_{xxz}$ tensor component of the quadratic ME response. This should not
be confused with the other second-order ME response, represented by the tensor
$\beta_{ijk}$ discussed in the previous section. In fact, $\beta_{ijk}$
identifies the magnetization induced at second order by an applied electric
field, whereas $\gamma_{ijk}$ corresponds to the polarization induced at
second order by an applied magnetic field. In the C-AFM phase of BiCoO3
$\beta_{ijk}$ is allowed only locally, i.e., separately in each sublattice,
but vanishes globally because of the time-reversal symmetry combined with
fractional translation. On the other hand, $\gamma_{ijk}$ is globally allowed
by symmetry because the lattice of BiCoO3 breaks inversion symmetry.

###  V.2 Ferromagnetic phase: giant magnetoelectric response

Figure 5: Plot of the change of polarization along the $x$ direction (total in
red circles and electronic response only in black squares) of the FM phase of
BiCoO3 under an applied magnetic field along the $x$ direction, hence probing
the $\alpha_{xx}$ ME tensor response.

As we have seen above, the total linear ME response of BiCoO3 is zero because
the local responses of each magnetic sublattice perfectly cancel each other
due to the C-AFM order. Here we explore by a computer experiment the ME
response of BiCoO3 with an artificial FM order. The FM order does not lead to
any crystal unit-cell doubling because it breaks also the time-reversal
symmetry, which (combined with a fractional translation) was instead preserved
in the C-AFM phase. Hence, we expect a net non-zero ME response in this case.
Assuming spins aligned along the z direction, the FM phase is described by the
magnetic space group $P4m^{\prime}m^{\prime}$, which allows for a diagonal ME
tensor with $\alpha_{xx}=\alpha_{yy}\neq\alpha_{zz}$.

In Fig. 5 we report our calculated ME response for the FM phase under an
applied Zeeman magnetic field along the $x$ direction (as obtained with the
ABINIT code) where the electronic response (clamped ions) and the full
response is shown. The full response is very large and gives a ME response
$\alpha_{xx}=1048$ ps/m, i.e., two orders of magnitude larger than what we
found for the individual ferromagnetic chains of the C-AFM order. Remarkably
this response is among the largest bulk (i.e., not composites [60]) crystal’s
ME response reported so far in type I multiferroics: among them we mention the
ME response of 480 ps/m for Fe2Mo3O8 [61], 730 ps/m for TbPO4 [3],
approximately 300 ps/m for DyFeO3 (beyond the critical magnetic field of 2.5T
at 3 K) [62], and approximately 1300 ps/m for Ni3TeO6 (at 2K) [63]. The
electronic contribution to the response is much smaller, about 16 ps/m, but it
is nevertheless one order of magnitude larger than those previously reported
in the literature: 0.34 ps/m for Cr2O3, 1.1 ps/m for LiNiPO4 [64], 1.8 ps/m
for FeS [65], and 1.8 ps/m for epitaxially strained CaMnO3 [20].

The large electronic response can be attributed to the fact that the band gap
is significantly reduced in the FM phase from 1.42 eV in the C-AFM phase to
0.82 eV, giving a large electronic contribution to the permittivity
($\varepsilon_{\infty}$) [20]. The large spin contribution to the response can
be attributed to the large magnetic susceptibility of the FM phase with
respect to the AFM phases (which increases the upper bound for the ME
response, Eq. (3)). Alternatively, it could be explained by the fact that when
applying an electric field, the only energy term to counteract in a FM phase
is the magnetocrystalline anisotropy energy (MAE), calculated to be
approximately 3 meV in our case, whereas in the AFM phase superexchange
interactions, which are much larger than the MAE (about 31 meV in the present
case), must be kept into account as well.

We cross-checked our result by computing the lattice-mediated ME response
using Eq. (6). To this end, we note that looking at the DMCs of the FM phase,
we find an enormous in-plane spin DMC value for all atoms, as follows: Bi atom
with 64 $\times 10^{-2}\mu_{\text{B}}/\text{\r{A}}$, Co atom with 1018 $\times
10^{-2}\mu_{\text{B}}/\text{\r{A}}$, apical O with $-439\times
10^{-2}\mu_{\text{B}}/\text{\r{A}}$, and basal O with 112 $\times
10^{-2}\mu_{\text{B}}/\text{\r{A}}$. The orbital contribution to the DMC,
although smaller, is significant for Co with 70 $\times
10^{-2}\mu_{\text{B}}/\text{\r{A}}$ and apical O with $-24\times
10^{-2}\mu_{\text{B}}/\text{\r{A}}$. Altogether, using the finite difference
technique, we obtain a lattice contribution to the ME response $\alpha_{xx}$
of 1115 ps/m, with a spin contribution of 1104 ps/m, which is in good
agreement with our applied magnetic field calculations. The out-of-plane
component, $\alpha_{zz}$, is found to be non-zero but considerably smaller,
around $-7$ ps/m.

The ground state of BiCoO3 is not FM, but, if one would be able to stabilize
BiCoO3 into a FM phase, it could potentially reveal a very large ME response.
In the conclusion, we discuss some strategies to do so.

##  VI Conclusions

We have studied the ME properties and the DMCs of bulk BiCoO3 in its $P4mm$
supertetragonal ground state structure. We found that the DMCs have both spin
and orbital origins, with the orbital contribution larger than the spin one in
the magnetic ground state; hence the orbital contribution cannot be neglected.
We also found that these DMCs are much larger than, e.g., the reference Cr2O3
ME crystal. We, however, established that BiCoO3 is not ME but rather anti-ME
as the ME responses of the two spin sublattices of the C-AFM ground state
phase are non-zero but opposite, resulting in a zero total ME response. We
calculated that each ME sublattice carries a rather large ME response
amplitude of about 85 ps/m and, as expected from the DMCs, the orbital
contribution is significantly larger (roughly 4 times larger) when compared to
the spin one.

Because each ferromagnetic sublattice of the C-AFM carries a large ME
response, we have computed this response in the higher energy FM phase of
BiCoO3 and found a giant response, of 1115 ps/m. We argue that the FM phase
has a much larger ME response than the single ME sublattice in the C-AFM phase
based on the fact that in the FM phase the applied magnetic field needs to
overcome the small MAE, whereas in the C-AFM phase an extra barrier is added
by the large AFM superexchange.

Our work indicates that by reaching the FM state in BiCoO3 a very large ME
response is unveiled. Considering this, let us discuss some of these
possibilities. One possibility is to uncompensate the C-AFM order through
doping, i.e., inducing ferrimagnetism, or to advance the material to a double
perovskite with two different B-sites. However, the superexchange barrier to
be overcome will remain the same in the ferrimagnetic system, although a large
ME response would still be possible, as shown by the response of each C-AFM
magnetic sublattice. Another problem would be to preserve the polar
distortions and the large magnetic exchange giving rise to room temperature
magnetoelectricity [66]. At last, doping may weaken the strong in-plane
coupling between Co atoms that occurs in the C-AFM phase, thus leading to a
decrease of $T_{\text{N}}$.

Another possibility is to grow BiCoO3 into superlattices or thin films [67],
aiming at making the FM phase energetically favorable, as in, for example,
LaMnO3 [68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78] or LaCoO3 [79]. A good
candidate for superlattices would be the Bi-based FM crystal BiMnO3 [80]. As a
possible route to switch to a FM order, we also mention the anionic
substitution, as it was done for, e.g., NaMnF3 [81], La0,7Ca0,3MnO3-δFx [82],
and SrNbO3-xNx [83]. We also stress that ferromagnetism can be produced at the
domain wall between AFM states [84], hence locally unveiling a large ME
response. Overall, we hope that our findings and suggestions motivate further
investigation of the possible large ME effect in BiCoO3 and, possibly, similar
compounds. To this end, we mention that a similar type of response could be
present in the twin system PbVO3 [85, 86, 87] as it shares the same
supertetragonal ferroelectric structure as BiCoO3 and the same magnetic order,
although with different chemistry and orbital filling of the transition metal
cation.

In this article, we have also provided an easy way for a 3D visualization of
the DMCs and to identify their multipolar nature. This is done by plotting, as
a vector field, the magnetic moment induced by displacing each atom
separately. The shape of the vector field provides useful insights on the
multipolar character of the DMCs. In this way we have been able to find that
the Bi DMC has a $x^{2}-y^{2}$ quadrupolar character, the Co and apical O DMCs
have both a monopolar and a $z^{2}$ quadrupolar origin, and the basal O DMC
has a toroidal and $xz$ and $yz$ quadrupolar nature. We confirmed these
results, which could be guessed from the aforementioned visual representation
of the DMCs, by directly computing from first principles the multipoles of the
atomic magnetization density.

Finally, we also emphasize that, in spite of the complexity of the
calculations (non-collinear magnetism within DFT+U with spin-orbit coupling,
which gives rise to multiple energy minima [64, 88, 89, 90], and the presence
of Co3+ cation in its high-spin state) and the small quantities at play (the
calculations of the DMCs from finite differences demand to probe magnetization
changes ranging from 10-5 to 10-2 $\mu_{\text{B}}$), we obtain very consistent
qualitative and quantitative results with two different codes (VASP and
ABINIT), i.e., with different implementations and different pseudopotentials,
in line with what was reported recently for the simpler case of Cr2O3 [91].
Remarkably, such an agreement is not initially guaranteed by using the codes
with a standard computational setup, i.e., with their usual default
parameters. Indeed, we obtain good convergence and agreements only if the
calculations are done with high care concerning not only the numerical
convergence parameters (plane wave cut-off, k-point grid for the reciprocal
space), but particularly the convergence of the wave functions during the
self-consistent cycle due to slow convergence of the magnetization. Our work,
therefore, shows that DFT codes are nowadays reliable enough and precise
enough to tackle such very subtle materials properties.

## Acknowledgements

We thank X. Gonze for fruitful discussions. E.B. and B.G. acknowledges the
FNRS and the Excellence of Science program (EOS “ShapeME”, No. 40007525)
funded by the FWO and F.R.S.-FNRS. EB and MB acknowledges the FNRS CDR project
“MULAN” No. (J.0020.20). Computational resources have been provided by the
Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de
la Recherche Scientifique (F.R.S.-FNRS) under Grant No. 2.5020.11 and the
Tier-1 Lucia supercomputer of the Walloon Region, infrastructure funded by the
Walloon Region under the grant agreement n°1910247 and by the High Performance
Computing Mesocenter of the University of Lille financed by the University,
the Hauts-de-France Region, the State, the FEDER and the University’s
laboratories through a pooling process.

## References

  * Fiebig [2005] M. Fiebig, Journal of Physics D: Applied Physics 38, R123 (2005).
  * Eerenstein  _et al._ [2006] W. Eerenstein, N. D. Mathur, and J. F. Scott, Nature 442, 759 (2006).
  * Rivera [2009] J.-P. Rivera, The European Physical Journal B 71, 299 (2009).
  * Kosub  _et al._ [2017] T. Kosub, M. Kopte, R. Hühne, P. Appel, B. Shields, P. Maletinsky, R. Hübner, M. O. Liedke, J. Fassbender, O. G. Schmidt, and D. Makarov, Nature Communications 8, 13985 (2017).
  * Kopyl  _et al._ [2021] S. Kopyl, R. Surmenev, M. Surmeneva, Y. Fetisov, and A. Kholkin, Materials Today Bio 12, 100149 (2021).
  * Liang  _et al._ [2021] X. Liang, H. Chen, and N. X. Sun, APL Materials 9, 041114 (2021).
  * Spaldin and Ramesh [2019] N. A. Spaldin and R. Ramesh, Nature Materials 18, 203 (2019).
  * Röntgen [1888] W. C. Röntgen, Annalen der Physik 271, 264 (1888).
  * Curie [1894] P. Curie, Journal de Physique Théorique et Appliquée 3, 393 (1894).
  * Debye [1926] P. Debye, Zeitschrift für Physik 36, 300 (1926).
  * Dzyaloshinskiǐ [1960] I. E. Dzyaloshinskiǐ, Soviet physics, JETP 10, 628 (1960).
  * Astrov [1960a] D. N. Astrov, Journal of Experimental and Theoretical Physics 11, 984 (1960a).
  * Astrov [1960b] D. N. Astrov, Journal of Experimental and Theoretical Physics 13, 1035 (1960b).
  * Folen  _et al._ [1961] V. J. Folen, G. T. Rado, and E. W. Stalder, Physical Review Letters 6, 607 (1961).
  * Rado and Folen [1961] G. T. Rado and V. J. Folen, Physical Review Letters 7, 310 (1961).
  * Brown  _et al._ [1968] W. F. Brown, R. M. Hornreich, and S. Shtrikman, Physical Review 168, 574 (1968).
  * Bayaraa  _et al._ [2021] T. Bayaraa, Y. Yang, M. Ye, and L. Bellaiche, Physical Review B 103, L060103 (2021).
  * Sasani  _et al._ [2022] A. Sasani, J. Íñiguez, and E. Bousquet, Physical Review B 105, 064414 (2022).
  * Fernández-Posada  _et al._ [2016] C. M. Fernández-Posada, A. Castro, J.-M. Kiat, F. Porcher, O. Peña, M. Algueró, and H. Amorín, Nature Communications 7, 12772 (2016).
  * Bousquet and Spaldin [2011] E. Bousquet and N. Spaldin, Physical Review Letters 107, 197603 (2011).
  * Diéguez and Íñiguez [2011] O. Diéguez and J. Íñiguez, Physical Review Letters 107, 057601 (2011).
  * Chun  _et al._ [2010] S. H. Chun, Y. S. Chai, Y. S. Oh, D. Jaiswal-Nagar, S. Y. Haam, I. Kim, B. Lee, D. H. Nam, K.-T. Ko, J.-H. Park, J.-H. Chung, and K. H. Kim, Physical Review Letters 104, 037204 (2010).
  * Malashevich  _et al._ [2012] A. Malashevich, S. Coh, I. Souza, and D. Vanderbilt, Physical Review B 86, 094430 (2012).
  * Íñiguez [2008] J. Íñiguez, Physical Review Letters 101, 117201 (2008).
  * Birol  _et al._ [2012] T. Birol, N. A. Benedek, H. Das, A. L. Wysocki, A. T. Mulder, B. M. Abbett, E. H. Smith, S. Ghosh, and C. J. Fennie, Current Opinion in Solid State and Materials Science 16, 227 (2012), multiferroics.
  * Bousquet and Cano [2016] E. Bousquet and A. Cano, Journal of Physics: Condensed Matter 28, 123001 (2016).
  * Fiebig  _et al._ [2016] M. Fiebig, T. Lottermoser, D. Meier, and M. Trassin, Nature Reviews Materials 1, 16046 (2016).
  * Xu  _et al._ [2024] C. Xu, H. Yu, J. Wang, and H. Xiang, Annual Review of Condensed Matter Physics 15, 85 (2024).
  * Ye and Vanderbilt [2014] M. Ye and D. Vanderbilt, Physical Review B 89, 064301 (2014).
  * Scaramucci  _et al._ [2012] A. Scaramucci, E. Bousquet, M. Fechner, M. Mostovoy, and N. A. Spaldin, Physical Review Letters 109, 197203 (2012).
  * Solovyev and Kolodiazhnyi [2016] I. V. Solovyev and T. V. Kolodiazhnyi, Physical Review B 94, 094427 (2016).
  * Yanagi  _et al._ [2018] Y. Yanagi, S. Hayami, and H. Kusunose, Physical Review B 97, 020404 (2018).
  * Ederer and Spaldin [2007] C. Ederer and N. A. Spaldin, Physical Review B 76, 214404 (2007).
  * Spaldin  _et al._ [2013] N. A. Spaldin, M. Fechner, E. Bousquet, A. Balatsky, and L. Nordström, Physical Review B 88, 094429 (2013).
  * Urru and Spaldin [2022] A. Urru and N. A. Spaldin, Annals of Physics 447, 168964 (2022).
  * Verbeek  _et al._ [2023] X. H. Verbeek, A. Urru, and N. A. Spaldin, Physical Review Research 5, L042018 (2023).
  * Note [1] Note that the phonon angular frequencies $\omega$ are related to the frequencies $\nu$, commonly expressed in THz or cm-1, by $\omega=2\pi\nu$.
  * Gonze and Lee [1997] X. Gonze and C. Lee, Physical Review B 55, 10355 (1997).
  * Gajdoš  _et al._ [2006] M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller, and F. Bechstedt, Physical Review B 73, 045112 (2006).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Computational Materials Science 6, 15 (1996).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Physical Review B 47, 558 (1993).
  * Gonze et al. [2020] X. Gonze et al., Computer Physics Communications 248, 107042 (2020).
  * Resta [2010] R. Resta, Journal of Physics: Condensed Matter 22, 123201 (2010).
  * Ceresoli  _et al._ [2010] D. Ceresoli, U. Gerstmann, A. P. Seitsonen, and F. Mauri, Phys. Rev. B 81, 060409 (2010).
  * THONHAUSER [2011] T. THONHAUSER, International Journal of Modern Physics B 25, 1429–1458 (2011).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Physical Review Letters 100, 136406 (2008).
  * Blöchl [1994] P. E. Blöchl, Physical Review B 50, 17953 (1994).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Physical Review B 13, 5188 (1976).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Physical Review B 23, 5048 (1981).
  * Liechtenstein  _et al._ [1995] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Physical Review B 52, R5467 (1995).
  * Ravindran  _et al._ [2008] P. Ravindran, R. Vidya, O. Eriksson, and H. Fjellvåg, Advanced Materials 20, 1353 (2008).
  * Uratani  _et al._ [2005] Y. Uratani, T. Shishidou, F. Ishii, and T. Oguchi, Japanese Journal of Applied Physics 44, 7130 (2005).
  * McLeod  _et al._ [2010] J. A. McLeod, Z. V. Pchelkina, L. D. Finkelstein, E. Z. Kurmaev, R. G. Wilks, A. Moewes, I. V. Solovyev, A. A. Belik, and E. Takayama-Muromachi, Physical Review B 81, 144103 (2010).
  * Bultmark  _et al._ [2009] F. Bultmark, F. Cricchio, O. Grånäs, and L. Nordström, Physical Review B 80, 035121 (2009).
  * Note [2] See https://github.com/materialstheory/multipyles.
  * Torrent  _et al._ [2008] M. Torrent, F. Jollet, F. Bottin, G. Zérah, and X. Gonze, Computational Materials Science 42, 337 (2008).
  * Jollet  _et al._ [2014] F. Jollet, M. Torrent, and N. Holzwarth, Computer Physics Communications 185, 1246 (2014).
  * Belik  _et al._ [2006] A. A. Belik, S. Iikubo, K. Kodama, N. Igawa, S.-i. Shamoto, S. Niitaka, M. Azuma, Y. Shimakawa, M. Takano, F. Izumi, and E. Takayama-Muromachi, Chemistry of Materials 18, 798 (2006).
  * Mostovoy  _et al._ [2010] M. Mostovoy, A. Scaramucci, N. A. Spaldin, and K. T. Delaney, Physical Review Letters 105, 087202 (2010).
  * Vaz  _et al._ [2010] C. A. F. Vaz, J. Hoffman, C. H. Ahn, and R. Ramesh, Advanced Materials 22, 2900 (2010).
  * Chang  _et al._ [2023] Y. Chang, Y. Weng, Y. Xie, B. You, J. Wang, L. Li, J.-M. Liu, S. Dong, and C. Lu, Physical Review Letters 131, 136701 (2023).
  * Tokunaga  _et al._ [2008] Y. Tokunaga, S. Iguchi, T. Arima, and Y. Tokura, Physical Review Letters 101, 097205 (2008).
  * Oh  _et al._ [2014] Y. S. Oh, S. Artyukhin, J. J. Yang, V. Zapf, J. W. Kim, D. Vanderbilt, and S.-W. Cheong, Nature Communications 5, 3201 (2014).
  * Bousquet and Spaldin [2010] E. Bousquet and N. Spaldin, Physical Review B 82, 220402 (2010).
  * Ricci and Bousquet [2016] F. Ricci and E. Bousquet, Physical Review Letters 116, 227601 (2016).
  * Kim  _et al._ [2019] M. K. Kim, J. Y. Moon, S. H. Oh, D. G. Oh, Y. J. Choi, and N. Lee, Scientific Reports 9, 5456 (2019).
  * Li  _et al._ [2022] P. Li, X.-S. Zhou, and Z.-X. Guo, npj Computational Materials 8, 20 (2022).
  * Gibert  _et al._ [2012] M. Gibert, P. Zubko, R. Scherwitzl, J. Íñiguez, and J.-M. Triscone, Nature Materials 11, 195 (2012).
  * Hou  _et al._ [2014] Y. S. Hou, H. J. Xiang, and X. G. Gong, Physical Review B 89, 064415 (2014).
  * Gibert  _et al._ [2015] M. Gibert, M. Viret, A. Torres-Pardo, C. Piamonteze, P. Zubko, N. Jaouen, J. M. Tonnerre, A. Mougin, J. Fowlie, S. Catalano, A. Gloter, O. Stéphan, and J. M. Triscone, Nano Letters 15, 7355 (2015).
  * An  _et al._ [2017] M. An, Y. Weng, H. Zhang, J.-J. Zhang, Y. Zhang, and S. Dong, Physical Review B 96, 235112 (2017).
  * Niu  _et al._ [2018] W. Niu, W. Liu, M. Gu, Y. Chen, X. Zhang, M. Zhang, Y. Chen, J. Wang, J. Du, F. Song, X. Pan, N. Pryds, X. Wang, P. Wang, Y. Xu, Y. Chen, and R. Zhang, Advanced Electronic Materials 4, 1800055 (2018).
  * Zhong  _et al._ [2018] C. Zhong, X. Lu, Y. Wan, Y. Min, Z. Zhao, P. Zhou, Z. Dong, and J. Liu, Journal of Magnetism and Magnetic Materials 466, 406 (2018).
  * Guan  _et al._ [2019] X. Guan, X. Shen, J. Zhang, W. Wang, J. Zhang, H. Wang, W. Wang, Y. Yao, J. Li, C. Gu, J. Sun, and R. Yu, Physical Review B 100, 014427 (2019).
  * Liu  _et al._ [2019] Y. Liu, H. Wong, K. Lam, C. Mak, and C. Leung, Journal of Magnetism and Magnetic Materials 481, 85 (2019).
  * Yao  _et al._ [2021] H. Yao, K. Jin, Z. Yang, Q. Zhang, W. Ren, S. Xu, M. Yang, L. Gu, E.-J. Guo, C. Ge, C. Wang, X. Xu, D. Zhang, and G. Yang, Advanced Materials Interfaces 8, 2101499 (2021).
  * Yang  _et al._ [2021] M. Yang, K. Jin, H. Yao, Q. Zhang, Y. Ji, L. Gu, W. Ren, J. Zhao, J. Wang, E.-J. Guo, C. Ge, C. Wang, X. Xu, Q. Wu, and G. Yang, Advanced Science 8, 2100177 (2021).
  * Lu  _et al._ [2022] Q. Lu, Z. Liu, Q. Yang, H. Cao, P. Balakrishnan, Q. Wang, L. Cheng, Y. Lu, J.-M. Zuo, H. Zhou, P. Quarterman, S. Muramoto, A. J. Grutter, H. Chen, and X. Zhai, ACS Nano 16, 7580 (2022).
  * Fuchs  _et al._ [2008] D. Fuchs, E. Arac, C. Pinta, S. Schuppler, R. Schneider, and H. v. Löhneysen, Physical Review B 77, 014434 (2008).
  * Choi  _et al._ [2017] E.-M. Choi, J. E. Kleibeuker, and J. L. MacManus-Driscoll, Scientific Reports 7, 43799 (2017).
  * Ren  _et al._ [2021] W.-n. Ren, K. Jin, E.-J. Guo, C. Ge, C. Wang, X. Xu, H. Yao, L. Jiang, and G. Yang, Physical Review B 104, 174428 (2021).
  * Altintas  _et al._ [2014] S. Altintas, N. Mahamdioua, A. Amira, and C. Terzioglu, Journal of Magnetism and Magnetic Materials 368, 111 (2014).
  * Garcia-Castro  _et al._ [2021] A. C. Garcia-Castro, Y. Ma, Z. Romestan, E. Bousquet, C. Cen, and A. H. Romero, Advanced Functional Materials 32, 10.1002/adfm.202107135 (2021).
  * Hirose  _et al._ [2017] H. T. Hirose, J.-i. Yamaura, and Z. Hiroi, Scientific Reports 7, 42440 (2017).
  * Ga  _et al._ [2024] Y. Ga, D. Yu, Y. Li, J. Jiang, L. Wang, P. Li, J. Liang, H. Feng, Y. Shi, and H. Yang, Physical Review B 109, L060402 (2024).
  * Bouvier  _et al._ [2023] P. Bouvier, A. Sasani, E. Bousquet, M. Guennou, and J. A. Moreira, Journal of Physics and Chemistry of Solids 173, 111092 (2023).
  * Shpanchenko  _et al._ [2004] R. V. Shpanchenko, V. V. Chernaya, A. A. Tsirlin, P. S. Chizhov, D. E. Sklovsky, E. V. Antipov, E. P. Khlybov, V. Pomjakushin, A. M. Balagurov, J. E. Medvedeva, E. E. Kaul, and C. Geibel, Chemistry of Materials 16, 3267 (2004).
  * Payne  _et al._ [2018] A. Payne, G. Avendaño-Franco, E. Bousquet, and A. H. Romero, Journal of Chemical Theory and Computation 14, 4455 (2018).
  * Payne  _et al._ [2019] A. Payne, G. Avedaño-Franco, X. He, E. Bousquet, and A. H. Romero, Physical Chemistry Chemical Physics 21, 21932 (2019).
  * Ponet  _et al._ [2023] L. Ponet, E. D. Lucente, N. Marzari, and et al., Research Square 10.21203/rs.3.rs-3358581/v1 (2023).
  * Bousquet  _et al._ [2024] E. Bousquet, E. Lelièvre-Berna, N. Qureshi, J.-R. Soh, N. A. Spaldin, A. Urru, X. H. Verbeek, and S. F. Weber, Journal of Physics: Condensed Matter 36, 155701 (2024).