# Characterizing and Overcoming Surface Paramagnetism in Magnetoelectric
Antiferromagnets

Sophie F. Weber  Materials Theory, ETH Zürich, Wolfgang-Pauli-Strasse 27, 8093
Zürich, Switzerland  Nicola A. Spaldin  Materials Theory, ETH Zürich,
Wolfgang-Pauli-Strasse 27, 8093 Zürich, Switzerland

###### Abstract

We use a combination of density functional theory and Monte Carlo calculations
to calculate the surface magnetization in magnetoelectric
$\mathrm{Cr_{2}O_{3}}$ at finite temperatures. Such antiferromagnets, lacking
both inversion and time-reversal symmetries, are required by symmetry to
posses an uncompensated magnetization density on particular surface
terminations. Here, we first show that the uppermost layer of magnetic moments
on the $(001)$ surface remain paramagnetic at the bulk Néel temperature,
bringing the theoretical estimate of surface magnetization density in line
with experiment. We demonstrate that the lower surface ordering temperature
compared to bulk is a generic feature of surface magnetization when the
termination reduces the effective Heisenberg coupling. We then propose two
methods by which the surface magnetization in $\mathrm{Cr_{2}O_{3}}$ could be
stabilised at higher temperatures. Specifically, we show that the effective
coupling of surface magnetic ions can be drastically increased either by a
different choice of surface Miller plane, or by $\mathrm{Fe}$ doping. Our
findings provide an improved understanding of surface magnetization properties
in AFMs.

Magnetoelectric (ME) antiferromagnets (AFMs) acquire a net magnetization
$\mathbf{M}$ in response to an applied electric field $\mathbf{E}$, and
conversely, a net electric polarization $\mathbf{P}$ in response to an applied
magnetic field $\mathbf{H}$Astrov (1960). For the linear ME effect to
manifest, an AFM must lack both inversion and time-reversal symmetries. This
symmetry requirement implies another intriguing property of ME AFMs; namely,
that certain surfaces must have a finite magnetic dipole densityBelashchenko
(2010). Such surface magnetization in ME AFMs has promising device
applications, since the ME effect allows the bulk domain to be readily
switched using electric fields in a constant magnetic fieldYe (2022), and the
direction of surface magnetization, which couples to the bulk AFM order
parameter, can be directly detectedHedrich  _et al._ (2021). Additionally,
surface magnetization plays a role in exchange bias coupling, extensively
exploited in magnetic sensors and storage devices to pin the magnetization
orientation of a ferromagnet (FM) by an adjacent AFMNogue´s and Schuller
(1999); Stamps (2000).
An important question about surface magnetism is its degree of disorder close
to the bulk Néel temperature $\mathrm{T}_{N}^{\mathrm{bulk}}$. Indeed, in the
case of $\mathrm{Cr_{2}O_{3}}$ (chromia), a prototypical ME AFM viewed as a
promising spintronics candidate due to its high Néel temperature of $\sim 300$
$\mathrm{K}$Echtenkamp (2021); Schlitz  _et al._ (2018); Muduli  _et al._
(2021); Ye (2022), theoretical predictions assuming that the bulk AFM order
persists at the surface greatly overestimate the size of the $(001)$ surface
magnetization density measured using nitrogen vacancy magnetometryAppel  _et
al._ (2019); Wörnle  _et al._ (2021); Spaldin (2021). This discrepancy is
resolved if the outermost surface $\mathrm{Cr}$ moments are disordered at the
measurement temperature, at or just below
$\mathrm{T}_{N}^{\mathrm{bulk}}$Spaldin (2021). In general, a better
understanding of the temperature dependence of surface magnetization in AFMs
would facilitate quantitative comparison between theory and experiment, and
could inform design of related spintronics devices.
In this letter, we use a combination of density functional theory (DFT) and
Monte Carlo (MC) calculations to explore the temperature dependence of surface
magnetism, taking $\mathrm{Cr_{2}O_{3}}$ as an example. We show that partial
to full disorder is a generic property of surface magnetization around the
bulk ordering temperature when surface magnetic moments have fewer or smaller
magnetic interactions than the bulk. We then propose two promising options for
stabilizing the surface magnetization of $\mathrm{Cr_{2}O_{3}}$ at the bulk
Néel temperature, first by using a Miller plane with a magnetic coupling close
to bulk, and secondly by adding a monolayer of $\mathrm{Fe}$ on the $(001)$
surface.
We first restate two key concepts, discussed in detail elsewhereStengel
(2011); Spaldin (2021). The first regards the construction of an
electrostatically stable, nonpolar surface termination for a given Miller
plane $(h,k,l)$Stengel (2011). A stable surface must have no bound charge,
since a finite $\sigma_{\mathrm{surf}}$ implies a diverging electrostatic
potentialNakagawa  _et al._ (2006). $\sigma_{\mathrm{surf}}$ is determined by
the component of bulk electric polarization $\mathbf{P}_{\mathrm{bulk}}$
perpendicular to the surfaceVanderbilt and King-Smith (1993):
$\mathbf{P}_{\mathrm{bulk}}\cdot\hat{\mathbf{n}}=\sigma_{\mathrm{surf}}$,
where $\hat{\mathbf{n}}$ is the unit surface normal. The periodicity of a bulk
crystal implies that $\mathbf{P}_{\mathrm{bulk}}$ is only defined modulo a
“polarization quantum” which corresponds to translating one electron by a
lattice vectorKing-Smith and Vanderbilt (1993). However, selecting a specific
surface termination dictates a particular basis choice for the bulk unit cell
(that which periodically tiles the semi-infinite solid containing the surface
of interest), and hence a single value of $\mathbf{P}_{\mathrm{bulk}}$. In
this case, $\sigma_{\mathrm{surf}}$ is single-valued, and a stable surface
plane has
$\mathbf{P}_{\mathrm{bulk}}\cdot\hat{\mathbf{n}}=\sigma_{\mathrm{surf}}=0$.
The second point relates to the connection between the bulk ME multipolization
tensor and the surface magnetizationSpaldin (2021). The multipolization tensor
is defined formally as $\mathcal{M}_{ij}=1/V\int
r_{i}\mathbf{\mu}_{j}(\mathbf{r})d^{3}\mathbf{r}$, where $r_{i}$ is the
$i\mathrm{th}$ cartesian component of position, $\mu_{j}(\mathbf{r})$ is the
$j\mathrm{th}$ component of magnetization density at position $\mathbf{r}$,
and $V$ is the unit cell volume. $\mathcal{M}$ describes first-order asymmetry
in $\mathbf{\mu}(\mathbf{r})$ beyond the magnetic dipoleSpaldin  _et al._
(2013). For materials in which $\mu(\mathbf{r})$ is localized around magnetic
ions, $\mathcal{M}_{ij}$ can be approximated by:

$\mathcal{M}_{ij}=\frac{1}{V}\sum_{\alpha}r^{\alpha}_{i}m^{\alpha}_{j},$ (1)

where the sum is over magnetic ions in the unit cell, and
$\mathbf{m}^{\alpha}$ is the local magnetic moment of atom $\alpha$.
The requirements for $\mathcal{M}$ to have nonzero components, that is, broken
inversion and time-reversal symmetries, are identical to those for a nonzero
linear ME response. Since a surface normal $\hat{\mathbf{n}}$ and an electric
field $\mathbf{E}$ are both polar vectors, introducing a surface reduces the
symmetry in the same way as applying an electric field in the bulk; therefore,
ME AFMs must have nonzero surface magnetizationSpaldin (2021); Belashchenko
(2010); He  _et al._ (2010). By analogy with the surface charge density
$\sigma_{\mathrm{surf}}$ resulting from the bulk polarization
$\mathbf{P}_{\mathrm{bulk}}$, the bulk multipolization tensor $\mathcal{M}$
gives rise to a surface magnetic dipole density $\mu_{\mathrm{surf}}$Spaldin
(2021), with the component $\mathcal{M}_{ij}$ giving the $\hat{i}$-oriented
magnetization density on a surface whose normal is parallel to $\hat{j}$. Like
$\mathbf{P}_{\mathrm{bulk}}$, the components of $\mathcal{M}$ are defined
modulo a multipolization “increment”, corresponding to moving a magnetic ion
by one lattice vector. But again, once a specific Miller plane and atomic
termination are selected, the origin of the bulk unit cell and thus the value
of multipolization are fixed. Therefore, the surface magnetization associated
with $\mathcal{M}$ is a single-valued quantity in the limit of bulk-like order
of the surface magnetic moments (see supplement for further discussion on the
connection between $\mathcal{M}$ and surface magnetization).

Figure 1: (a) Hexagonal unit cell of $\mathrm{Cr_{2}O_{3}}$ with nearest neighbor Heisenberg couplings indicated for bulk (solid lines) and $(001)$ surface (dashed lines) $\mathrm{Cr}$s. Numbers in parentheses indicate the coupling degeneracies (not all nearest neighbors are visible in the diagram). Lighter colored atoms are in the adjacent unit cell. White numbers label the $12$ magnetic $\mathrm{Cr}$ in the $(001)$ slab for ease of discussion. (b) Left: unit cell which defines the nonpolar $(001)$ surface with ground state AFM magnetic ordering. Right: unit cell used to calculate surface magnetization if the top moment of the nonpolar surface is disordered. Table 1: Heisenberg coupling constants, with degeneracies in parentheses, and total effective coupling $\lambda_{J}^{i}$ calculated in this work for bulk $\mathrm{Cr_{2}O_{3}}$, for a $\mathrm{Cr}$ ion on a $(001)$ surface (using the bulk values and those computed from a relaxed slab respectively), for a $\mathrm{Cr}$ on a $(\bar{1}02)$ surface, and $\mathrm{Cr}$-$\mathrm{Fe}$ couplings for $[001]$ $\mathrm{Cr_{2}O_{3}}$ with an $\mathrm{Fe}$ monolayer. | bulk |  $(001)$ surface |  $(\bar{1}02)$ surface |  $\mathrm{Fe}$ monolayer on $(001)$ surface
---|---|---|---|---
|  | bulk | relaxed | bulk | relaxed | relaxed
$J_{1}$ ($\mathrm{meV}$) | -10.46 (1) | - | - | -10.46 (1) | -17.16 (1) | -
$J_{2}$ ($\mathrm{meV}$) | -7.88 (3) | - | - | -7.88 (2) | -9.42 (2) | -
$J_{3}$ ($\mathrm{meV}$) | +0.86 (3) | +0.86 (3) | -0.15 (3) | +0.86 (1) | +0.43 (1) | -30.81 (3)
$J_{4}$ ($\mathrm{meV}$) | +1.22 (6) | +1.22 (3) | +4.44 (3) | +1.22 (5) | +0.43a (4)/+3.40b (1) | -19.77 (3)
$J_{5}$ ($\mathrm{meV}$) | -1.41 (1) | -1.41 (1) | -0.39 (1) | - | - | -3.63 (1)
$\lambda_{J}$ ($\mathrm{meV}$) | 40.25 | 2.48 | 14.18 | 31.47 | 40.16 | 155.33

_Results and Discussion._ –$\mathrm{Cr_{2}O_{3}}$ crystallizes in the corundum
structure with magnetic space group $\mathrm{R\bar{3}^{\prime}c^{\prime}}$
[161]Fechner  _et al._ (2018). Figure 1(a) shows the $12$-$\mathrm{Cr}$ unit
cell in the hexagonal setting. Bulk $\mathrm{Cr_{2}O_{3}}$ adapts an “up down
up down” ordering of the $\mathrm{Cr}$ magnetic moments along $[001]$ as shown
for the unit cells in Figure 1(b). This ground state order is well
describedSamuelsent  _et al._ (1970); Shi  _et al._ (2009) by a Heisenberg
Hamiltonian,

$\mathcal{H}_{\mathrm{Heis}}=-\sum\limits_{<i,j>}J_{i,j}(\mathbf{e}_{i}\cdot\mathbf{e}_{j}),$
(2)

that includes coupling up to the fifth nearest neighbors, where
$\mathbf{e}_{i}$ is the unit vector parallel to the local magnetic moment of
the $\mathrm{Cr}$ ion at site $i$, and $J_{i,j}$ is the Heisenberg coupling
constant between spins $i$ and $j$. The couplings $J_{1}$-$J_{5}$, where
$J_{n}$ denotes the coupling for the $\mathrm{n}^{th}$ nearest neighbor, are
depicted in Figure 1(a). The quantitative values of $J_{1}$-$J_{5}$ for bulk
$\mathrm{Cr_{2}O_{3}}$, which we calculate with the method outlined in
reference 22 using first-principles DFT+U as implemented in the VASP
softwareKresse and Furthmüller (1996), are given in Table 1. Our values are in
good agreement with previous DFT calculations using similar parametersShi  _et
al._ (2009). The magnetism is dominated by the strong AFM $J_{1}$ and $J_{2}$
couplings.
We first review magnetism on the $(001)$ surface of vacuum-terminated chromia.
The bulk unit cell with a single terminating $\mathrm{Cr}$ on the left-hand
side of Figure 1(b) defines the nonpolar $(001)$ surface according to the
formula $\mathbf{P}_{\mathrm{bulk}}=1/V\sum_{i}Z_{i}\mathbf{r}_{i}$, where
$Z_{i}$ is the formal ionic charge ($+3$ and $-2$ for $\mathrm{Cr}$ and
$\mathrm{O}$ respectively), and $\mathbf{r}_{i}$ the position of atom $i$ in
the unit cell. If we assume all $\mathrm{Cr}$ magnetic moments are fully
polarized along $[001]$ with the bulk AFM order, using the formal value
$3\mathrm{\mu_{B}}$ for $\mathrm{Cr^{3+}}$ and the fractional coordinates
(given in the supplement) in the hexagonal cell, equation 1 yields a
$+\hat{z}/[001]$-oriented $(001)$ surface magnetization of $11.9$
$\mathrm{\mu_{B}}/\mathrm{nm^{2}}$ for the magnetic domain depicted (all other
components of the multipolization tensor are zero within the local moment
approximation; small $(1,1)$ and $(2,2)$ components are symmetry-allowed if
one uses the exact integral formUrru and Spaldin (2022)). The energetically
equivalent AFM domain in which the directions of all magnetic moments are
reversed has a value of equal magnitude and opposite sign.
As mentioned previously, this theoretical predication overestimates
measurements of $(001)$ $\mathrm{Cr_{2}O_{3}}$ surface magnetism using
scanning nitrogen vacancy magnetometryAppel  _et al._ (2019); Wörnle  _et al._
(2021), which yield values between $1.6$ to $2.3$
$\mathrm{\mu_{B}}/\mathrm{nm^{2}}$ (the sign of magnetization cannot be
directly determined). Recall however that the $11.9$
$\mathrm{\mu_{B}}/\mathrm{nm^{2}}$ value is calculated assuming that all
$\mathrm{Cr}$ magnetic moments are fully ordered along $[001]$. Looking at the
outermost $\mathrm{Cr}$ for the nonpolar termination in Figure 1(a), we see
that it lacks $J_{1}$ and $J_{2}$ nearest neighbors, only retaining the
smaller $J_{3}$-$J_{5}$ couplings. From a mean-field argument, the ordering
temperature for a given magnetic moment at site $i$ is proportional to
$\lambda^{i}_{J}S_{i}$Mostovoy  _et al._ (2010), where $S_{i}$ is the spin
value and the total effective Heisenberg coupling for site $i$ is

$\lambda^{i}_{J}=\sum\limits_{j}J_{ij}\times(\hat{e}_{i}\cdot\hat{e_{j}}).$
(3)

Using the $J$ values in Table 1 calculated for bulk $\mathrm{Cr_{2}O_{3}}$,
$\lambda_{J}$ for a bulk $\mathrm{Cr}$ spin is
$\lambda_{J}^{\mathrm{bulk}}=-J_{1}-3J_{2}-3J_{3}+6J_{4}-J_{5}=40.25\mathrm{meV}$,
whereas the $\mathrm{Cr}$ on the $(001)$ surface ($\mathrm{Cr}$ $12$ with the
convention in Figure 1(a)) has
$\lambda_{J}^{\mathrm{surf}}=-3J_{3}+3J_{4}-J_{5}=2.48\mathrm{meV}$. Thus,
$T_{N}^{\mathrm{surf}}/T_{N}^{\mathrm{bulk}}=0.06$, implying that for the room
temperature $T\sim 293\mathrm{K}$ magnetometry measurements, just below
$\mathrm{T_{N}^{bulk}}\sim 300\mathrm{K}$, we expect the surface $\mathrm{Cr}$
to be paramagnetic. Taking into account this magnetic dead layer, a more
appropriate basis for predicting the surface magnetization is that shown on
the right-hand side of Figure 1(b), corresponding to removing the surface
magnetic moment by displacing it downwards one $\mathbf{c}$ lattice vector.
Recalculating the $(3,3)$ component of $\mathcal{M}$ using the $\mathrm{Cr}$
positions of this new unit cell yields $\mu_{(001)}^{z}=-2.4$
$\mathrm{\mu_{B}}/\mathrm{nm^{2}}$, in good agreement with experiment.

Figure 2: (a) Normalized magnetization along $\hat{z}/\mathrm{[001]}$ as a
function of temperature for a bulk $\mathrm{Cr}$ in the center of a
$[001]$-oriented slab, and for the surface $\mathrm{Cr}$ (labeled on the
structure in the inset). Here the Heisenberg couplings calculated from bulk
$\mathrm{Cr_{2}O_{3}}$ (first column of Table 1) are used for both surface and
bulk $\mathrm{Cr}$ moments. (b) $m_{z}$ using couplings calculated from a
relaxed slab for the surface $\mathrm{Cr}$.

To confirm our analysis, we examine the temperature dependence of
magnetization by performing Monte Carlo (MC) simulations as implemented in the
UppASD spin dynamics packageSkubic  _et al._ (2008) of a $[001]$-oriented
$\mathrm{Cr_{2}O_{3}}$ slab using a $42\times 42\times 1$ supercell of the
$12$-$\mathrm{Cr}$-atom hexagonal unit cell, having checked that this
thickness, with six $\mathrm{Cr_{2}O_{3}}$ layers, is sufficient to capture
both bulk and surface behavior. We enforce in-plane periodic boundary
conditions and vacuum boundary conditions along $[001]$. Figure 2(a) shows the
absolute value of the $\hat{z}$ component of bulk magnetization as a function
of temperature, calculated by averaging the projected $m_{z}$ of the
$6^{\mathrm{th}}$ $\mathrm{Cr}$ sublattice in the center of our unit-cell
thick slab, compared to the averaged $m_{z}$ of the terminating $\mathrm{Cr}$s
on the nonpolar $(001)$ surface. We also confirmed that all sublattices other
than the outermost $\mathrm{Cr}$ on both sides of the slab have bulk-like
$m_{z}(T)$ behavior. We see that, whereas the center “bulk” $\mathrm{Cr}$
exhibits the normal Langevin-like $m_{z}(T)$ curve, the surface magnetization
falls off rapidly with increasing temperature and is negligible at
$T_{N}^{\mathrm{bulk}}$, consistent with earlier combined DFT-MC
calculationsWysocki  _et al._ (2012). Note that our calculated Heisenberg
constants lead to a significant underestimate of $T_{N}^{\mathrm{bulk}}$
($T_{N}^{\mathrm{bulk}}\sim 170$ $\mathrm{K}$ based on Figure 2); this has
been observed in previous DFT-MC calculations of $\mathrm{Cr_{2}O_{3}}$Kota
_et al._ (2013).
While $m_{z}(T)$ in Figure 2(a) for both surface and bulk $\mathrm{Cr}$s are
computed using the DFT $J$ values calculated with bulk $\mathrm{Cr_{2}O_{3}}$,
atomic relaxation can lead to significant renormalization of the surface
couplings. The third column of Table 1 shows the values of $J_{3}$, $J_{4}$
and $J_{5}$ for the surface $\mathrm{Cr}$ computed using a $[001]$ vacuum-
terminated $12$-$\mathrm{Cr}$-thick slab which we structurally relax within
DFT. The effective coupling for the surface $\mathrm{Cr}$ when taking
relaxation into account is
$\lambda_{J}^{\mathrm{surf,relaxed}}=14.2\mathrm{meV}$. Figure 2(b) shows
$m_{z}(T)$ for the surface $\mathrm{Cr}$ with these relaxed values (we keep
the remaining $J$s for the the ten non-surface $\mathrm{Cr}$ set to bulk
values, having checked that the coupling renormalization for these ions upon
relaxation negligibly affects the results). While the surface $m_{z}$ is still
disordered at $T_{N}^{\mathrm{bulk}}$, the increased
$\lambda_{J}^{\mathrm{surf}}$ leads to a roughly linear decrease of $m_{z}$
with increasing $T$, as opposed to the exponential-like falloff in Figure
2(a).

Figure 3:  Ratio of $m_{z}^{\mathrm{surf}}$ to $m_{z}^{\mathrm{bulk}}$ at
$T=100$ $\mathrm{K}$ ($\sim T_{N}^{\mathrm{bulk}}/2$) as a function of
$\lambda_{J}^{\mathrm{surf}}$. $\lambda_{J}^{\mathrm{surf}}$ is fixed in the
MC simulations by setting the three $J_{4}^{\mathrm{surf}}$ to
$\lambda_{J}^{\mathrm{surf}}/3$ and all other surface couplings to zero. The
dashed blue line shows the value of $\lambda_{J}^{\mathrm{bulk}}$.

To determine the dependence of the surface magnetism on
$\lambda_{J}^{\mathrm{surf}}$ in detail, we next vary
$\lambda_{J}^{\mathrm{surf}}$ manually in the MC simulations by fixing each of
the three $J_{4}^{\mathrm{surf}}$ to one third the desired
$\lambda_{J}^{\mathrm{surf}}$, while setting all other surface couplings to
zero. For each value of $\lambda_{J}^{\mathrm{surf}}$ we calculate
$\frac{m_{z}^{\mathrm{surf}}}{m_{z}^{\mathrm{bulk}}}$ at $T=100$ $\mathrm{K}$;
the result is plotted in Figure 3. We choose $100$ $\mathrm{K}$ as a
representative temperature because it is roughly $T_{N}^{\mathrm{bulk}}/2$
($\lambda_{J}^{\mathrm{surf}}$ marginally affects $T_{N}^{\mathrm{bulk}}$ due
to the finite slab size, thus $T_{N}^{\mathrm{bulk}}/2$ ranges from $85$ to
$110$ $\mathrm{K}$ for the range of $\lambda_{J}^{\mathrm{surf}}$ in Figure
3). $\frac{m_{z}^{\mathrm{surf}}}{m_{z}^{\mathrm{bulk}}}|_{T=100\mathrm{K}}$
increases monotonically with $\lambda_{J}^{\mathrm{surf}}$ and matches the
bulk magnetization ($\frac{m_{z}^{\mathrm{surf}}}{m_{z}^{\mathrm{bulk}}}=1$)
roughly when $\lambda_{J}^{\mathrm{surf}}$ equals $\lambda_{J}$ for the bulk
$\mathrm{Crs}$ (dashed blue line). Therefore, engineering
$\lambda_{J}^{\mathrm{surf}}$ to be close to $\lambda_{J}^{\mathrm{bulk}}$ can
be taken as a criterion for obtaining bulk-like temperature dependence of
surface magnetization. Moreover, if the Heisenberg $J$s for a material are
known, one can quickly calculate $\lambda_{J}^{\mathrm{surf}}$ and estimate
how much $m^{\mathrm{surf}}$ is likely to be reduced relative to bulk
magnetization.

Figure 4: (a) Bulk unit cell defining nonpolar $(\bar{1}02)$ surface of
$\mathrm{Cr_{2}O_{3}}$, with $[001]$ oriented bulk magnetic order. (b)
Heisenberg couplings for $\mathrm{Cr}$ ions on the $(\bar{1}02)$ surface. (c)
$|m_{z}(T)|$ (with $z$ along the bulk $[001]$ easy axis) for $\mathrm{Cr}$ in
bulk and on a $(\bar{1}02)$ surface (using relaxed surface couplings). (d)
$|m_{z}(T)|$ for $\mathrm{Fe}$, $\mathrm{Cr}$ $10$, and a central bulk
$\mathrm{Cr}$ (positions $12$, $10$ and $6$ respectively from Figure 1(a)).

We now discuss two approaches, which can also be applied to other ME AFMs, for
stabilizing surface magnetization in $\mathrm{Cr_{2}O_{3}}$ at higher
temperatures. The first is to use a surface corresponding to a different
Miller plane for which $\lambda_{J}^{\mathrm{surf}}$ is close to
$\lambda_{J}^{\mathrm{bulk}}$. We demonstrate this for $\mathrm{Cr_{2}O_{3}}$
in Figure 4. Figure 4(a) shows the unit cell corresponding to the non-polar
termination of a $(\bar{1}02)$ surface, which has a large
$\lambda_{J}^{\mathrm{surf}}$ and a non-negligible theoretical surface
magnetization density. Specifically, for the domain shown, Equation 1 predicts
an out-of-plane (in-plane) magnetization component of $-4.75$
$\mathrm{\mu_{B}}/\mathrm{nm^{2}}$ ($+7.55$
$\mathrm{\mu_{B}}/\mathrm{nm^{2}}$) respectively on the $(\bar{1}02)$ surface
(see supplement).
The couplings and degeneracies, shown in Figure 4(b), retained by the
outermost $\mathrm{Cr}$ and corresponding $\lambda_{J}^{\mathrm{surf}}$ for
the $(\bar{1}02)$ surface are given in Table 1 both with DFT values calculated
from bulk $\mathrm{Cr_{2}O_{3}}$ and with surface couplings calculated from a
relaxed $[\bar{1}02]$-oriented slab. The $a$ and $b$ superscripts refer to the
two $J_{4}$ couplings depicted in Figure 4(b) which become inequivalent upon
relaxation (we double the unit cell in the surface plane in order to show all
couplings). The overall $\lambda_{J}^{\mathrm{surf}}$ is $31.47$ ($40.16$)
$\mathrm{meV}$ for bulk (relaxed) coupling values. Even using the bulk values
for $(\bar{1}02)$ surface moments, the surface magnetization is nearly bulk-
like, and with the relaxed values leading to
$\lambda_{J}^{\mathrm{surf}}\sim\lambda_{J}^{\mathrm{bulk}}$, bulk and
$(\bar{1}02)$ surface $m_{z}(T)$ lie on top of each other (Figure 4(c)).
Extensive research has been devoted to the application of $[001]$
$\mathrm{Cr_{2}O_{3}}$ films in spintronic memory devices, where the AFM bulk
domain serves as a logical bit whose direction can be read out by the sign of
the surface magnetization (this is usually determined indirectly via the sign
of the hysteresis loop shift, i.e. exchange bias, in an adjacent FMBorisov
_et al._ (2005); He  _et al._ (2010); Ye (2022)). Our results imply that
magnetism on the $(\bar{1}02)$ surface in chromia is strongly coupled to the
underlying bulk AFM domain, even at $T_{N}^{\mathrm{bulk}}$, in contrast to
the $(001)$ surface where the surface is essentially paramagnetic at room
temperature. Thus, a $\mathrm{Cr_{2}O_{3}}$-based device with a $(\bar{1}02)$
rather than $(001)$ surface plane might be a more robust option for memory
applications. More fundamentally, a comparison of exchange bias properties for
the $(001)$ and $(\bar{1}02)$ surfaces could shed light on the underlying
mechanism.
Our second proposed method for stabilizing surface magnetization involves
chemical substitution. We take the $\mathrm{Cr_{2}O_{3}}$ $(001)$ surface, and
deposit a monolayer of $\mathrm{Fe}$ on top (substituting $\mathrm{Cr}$ in the
$12^{\mathrm{th}}$ position in Figure 1(a)). Since $\mathrm{Fe}$ adopts a $3+$
valence state, this structure is nonpolar and stable (see supplement for
further discussion). Crucially, while the $\mathrm{Cr}$-$\mathrm{Cr}$
$J_{3}$-$J_{5}$ are negligible compared to $J_{1}$ and $J_{2}$, prior DFT
studies using a $\mathrm{Cr_{2}O_{3}}$-$\mathrm{Fe_{2}O_{3}}$ heterostructure
indicate that the $J_{3}$ and $J_{4}$ $\mathrm{Cr}$-$\mathrm{Fe}$ couplings at
the interface are tens of $\mathrm{meV}$Kota  _et al._ (2014). The difference
in strengths and signs of $\mathrm{Cr}$-$\mathrm{Cr}$ and
$\mathrm{Cr}$-$\mathrm{Fe}$ couplings in the corundum structure can be
attributed to the relative $e_{g}$-$t_{2g}$ occupation of the
$\mathrm{Cr^{3+}}$ $(t_{2g}^{3},e_{g}^{0})$ and $\mathrm{Fe^{3+}}$
$(t_{2g}^{3},e_{g}^{2})$ ions, combined with the coupling angles via
oxygenKota  _et al._ (2013) In the final column of Table 1, we show our
results for the surface $\mathrm{Fe}$-$\mathrm{Cr}$ couplings calculated using
a relaxed $\mathrm{Cr_{2}O_{3}}$ slab terminated on one side with
$\mathrm{Fe}$. $J_{3}^{\mathrm{surf,Cr-Fe}}$ and $J_{4}^{\mathrm{surf,{}Cr-
Fe}}$ are even larger than the $J_{1}^{\mathrm{Cr-Cr}}$ and
$J_{2}^{\mathrm{Cr-Cr}}$ that are dominant in $\mathrm{Cr_{2}O_{3}}$ bulk.
Figure 4(d) shows the absolute value of $m_{z}(T)$ for the surface
$\mathrm{Fe}$ monolayer and center $\mathrm{Cr}$ bulk, as well as for
“$\mathrm{Cr}$ $10$” (according to the labeling in Figure 1(a)) which is
coupled to the $\mathrm{Fe}$ via $J_{4}^{\mathrm{surf,{}Fe-Cr}}$ (note that
$\mathrm{Cr}$ $10$ reverses its orientation from that in bulk due to the
strong AFM $J_{4}^{\mathrm{surf,{}Fe-Cr}}$ coupling). The $\mathrm{Cr}$s
coupled directly to the $\mathrm{Fe}$ monolayer have magnetization
intermediate between those of the deeper bulk $\mathrm{Cr}$ and of
$\mathrm{Fe}$.
A notable feature of the surface $\mathrm{Fe}$ monolayer magnetization in
Figure 4(d) is that $m_{z}^{\mathrm{Fe,surf}}$ (and the $\mathrm{Cr}$ directly
below) order at higher temperatures than the bulk $\mathrm{Cr}$, making it an
attractive test case for fundamental research in paramagnetic bulk materials
with surface magnetic orderDedkov  _et al._ (2007); Rosenberg and Franz
(2012); Schulz  _et al._ (2019). Moreover, if scanning nitrogen vacancy
magnetometry measurements of $\mathrm{Fe}$-capped $\mathrm{Cr_{2}O_{3}}$ could
be compared at temperatures just above (where only $\mathrm{Fe}$ and the top-
most $\mathrm{Cr}$ are ordered) and below $T_{N}^{\mathrm{bulk}}$, monitoring
how the measured surface magnetization changes would provide a clear
indication of the technique’s depth resolution.
In summary, we have examined finite temperature properties of surface
magnetization in AFMs using ME $\mathrm{Cr_{2}O_{3}}$ as an example. Our
combined DFT-MC calculations demonstrate that disorder of surface magnetic
ions at $T_{N}^{\mathrm{bulk}}$ likely explains the discrepancy between
theoretical and experimental surface magnetization estimates on $(001)$
$\mathrm{Cr_{2}O_{3}}$. We establish a framework for assessing the relative
ordering temperature of surface and bulk magnetization based on effective
Heisenberg couplings. Finally, we have discussed two options for stabilizing
surface magnetism, which would allow for higher temperature operation of
relevant spintronic devices. We hope this work stimulates efforts, both
theoretical and experimental, to better understand and characterize surface
magnetization in AFMs.

###### Acknowledgements.

We thank Xanthe Verbeek, Tara Tošić, Kai Wagner, Paul Lehman, Patrick
Maletinsky, Sayantika Bhowal and Andrea Urru for useful discussions. This work
was funded by the ERC under the European Union’s Horizon 2020 research and
innovation program project HERO with grant number 810451. Calculations were
performed at the Swiss National Supercomputing Centre (CSCS) under project
number s889 and on the EULER cluster of ETH Zürich.

##  Appendix A Density Functional Theory calculation details

In order to calculate the relaxed structures and Heisenberg coupling constants
for $\mathrm{Cr_{2}O_{3}}$ in this work, we use density functional theory
(DFT), employing the Vienna _ab initio_ simulation package (VASP)Kresse and
Furthmüller (1996) with the localized density approximation (LDA) within the
projector augmented wave method (PAW)Blochl (1994). We use the standard VASP
PAW pseudopotentials with the following valence electron configurations:
$\mathrm{Cr}$ $(3p^{6}4s^{2}3d^{5})$, $\mathrm{O}$ $(2s^{2}3p^{4})$, and
$\mathrm{Fe}$ $(4s^{2}3d^{6})$ (for calculations of $\mathrm{Fe}$-capped
$\mathrm{Cr_{2}O_{3}}$). We use collinear spin-polarized calculations,
neglecting spin-orbit coupling except when calculating the magnetocrystalline
anisotropy. We use an energy cutoff of $800$ $\mathrm{eV}$ for our plane wave
basis set, and a Gamma-centered $9\times 9\times 5$ k-mesh for the bulk
$30$-atom hexagonal unit cell of $\mathrm{Cr_{2}O_{3}}$. To model the $(001)$
($(\bar{1}02)$) $\mathrm{Cr_{2}O_{3}}$ surfaces we use hexagonal (monoclinic)
cells with $9\times 9\times 1$ ($11\times 6\times 1$) Gamma-centered k-meshes
with $15$ $\mathrm{\AA}$ vacuum in the direction of the surface normals. We
use the tetrahedron methodBlochl (1994) for Brillouin zone integrations. We
find that these parameters lead to total energy convergence of $<1$
$\mathrm{meV}$ per formula unit. We relax all structures, both bulk and slabs,
until forces on all atoms are less than $0.01$ $\mathrm{eV}/\mathrm{\AA}$.
To approximately capture the localized nature of $3d$ electrons in
$\mathrm{Cr}$, we add a Hubbard U correctionAnisimov  _et al._ (1997) using
the rotationally invariant method by Dudarev et al.Dudarev  _et al._ (1998).
We set $\mathrm{U}=4$ $\mathrm{eV}$ based on prior DFT+U work on
$\mathrm{Cr_{2}O_{3}}$ with $\mathrm{U}=4$, $\mathrm{J}\sim 0.5$
$\mathrm{eV}$Shi  _et al._ (2009); Fechner  _et al._ (2018) (we also use
$\mathrm{U}=4$ on the $\mathrm{Fe}$ $\mathrm{d}$ states for calculations of
the $\mathrm{Fe}$ capped $[001]$ chromia slab). We find that including the
Hund’s coupling $\mathrm{J}=0.5$ $\mathrm{eV}$ does not significantly affect
the computed Heisenberg coupling constants, thus we only use the Hubbard
$\mathrm{U}$. As mentioned in the main text, with this $\mathrm{U}$ value our
calculated couplings lead to an underestimated bulk Néel temperature
$T_{N}^{\mathrm{bulk}}$ from the Monte Carlo simulations, which was also the
case for a prior DFT-MC study of chromia using a similar U valueKota  _et al._
(2014). Using a smaller $\mathrm{U}$ ($\mathrm{U}=2$ $\mathrm{eV}$ for
reference 25) would lead to larger Heisenberg couplings (due to decreased
localisation), and hence a $T_{N}^{\mathrm{bulk}}$ closer to
experimentMostovoy  _et al._ (2010). However, we choose $\mathrm{U}=4$
$\mathrm{eV}$ because with this value we achieve the correct sign of
magnetocystalline anisotropy energy (MAE), i.e. easy axis along the hexagonal
$[001]$ directionDudko  _et al._ (1971); Tobia  _et al._ (2010). With
$\mathrm{U}=2$ $\mathrm{eV}$ on the other hand, we calculate a qualitatively
incorrect easy plane. Thus, we believe the higher $\mathrm{U}$ value overall
better describes the magnetic properties of $\mathrm{Cr_{2}O_{3}}$. We also
note that the _relative_ values of $J_{1}$ through $J_{5}$ are similar for a
wide range of $\mathrm{U}$ as can be seen from reference 21. Thus, the
primarily qualitative conclusions drawn in our work hold in spite of the
$T_{N}^{\mathrm{bulk}}$ underestimate.
We compute the Heisenberg couplings for both bulk and slab structures using
the method outlined in reference 22. Essentially, the coupling between two
specific sites $i$ and $j$ is calculated from four total energy calculations
in which the magnetic moments on these two sites are set to
$(i,j)=(\uparrow\uparrow),(\uparrow\downarrow),(\downarrow\uparrow),(\downarrow\downarrow)$,
with moments on all other sites in the unit cell kept constant. In the case of
the relaxed vacuum-terminated structures, this method allows us to calculate
the Heisenberg couplings for each $\mathrm{Cr}$ site and thus differentiate
between the values for the surface $\mathrm{Cr}$ and for $\mathrm{Cr}$ in the
center of the slab which retain couplings close to the results from bulk.

##  Appendix B Monte Carlo calculation details

To explore the temperature dependence of surface magnetism in
$\mathrm{Cr_{2}O_{3}}$, we use Monte Carlo (MC) simulations as implemented in
the UppASDSkubic  _et al._ (2008) spin dynamics package. We use supercells
with $42\times 42\times 1$ ($38\times 19\times 1$) magnetic unit cells for the
$(001)$ ($(\bar{1}02)$) surfaces, for a total of $21168$ ($23104$) magnetic
atoms in the simulation box. We use periodic boundary conditions in the in-
plane $\mathbf{a}$ and $\mathbf{b}$ directions and vacuum boundary conditions
in the out-of-plane $\mathbf{c}$ direction. To test convergence of our
results, we also performed MC simulations with a simulation box doubled along
the $\mathrm{[}001]$ surface normal, i.e. $24$-$\mathrm{Cr}$ tall as opposed
to the $12$-$\mathrm{Cr}$ tall box used in the main manuscript. We found that
the projected $m_{z}(T)$ for both surface $\mathrm{Cr}$ and $\mathrm{Cr}$
sublattices in the center of the slab did not change noticeably upon doubling
the height; thus, the $12$-$\mathrm{Cr}$ tall box ($16$-$\mathrm{Cr}$ tall for
the $(\bar{1}02)$ surface) is sufficiently thick to capture behavior of both
bulk and surface magnetization. To prevent the magnetization axis from
drifting, we add a uniaxial magnetoanisotropy energy (MAE) of $0.11$
$\mathrm{meV}$ per $\mathrm{Cr}$ along the $\hat{z}$/$[001]$ direction, along
the lines of previous studiesMostovoy  _et al._ (2010). The experimental MAE
of bulk $\mathrm{Cr_{2}O_{3}}$, as well as the value we calculate with DFT+U
(about $4$ $\mu\mathrm{eV}$ per $\mathrm{Cr}$) is two orders of magnitude
smallerDudko  _et al._ (1971); Tobia  _et al._ (2010). However, due to the
finite size of the simulation box in the MC simulations the MAE must be scaled
up to prevent unphysical fluctuations of the magnetization axis. Simulations
at each temperature were performed with $2\times 10^{5}$ initial steps to
bring the system to thermal equilibrium, and $N_{MC}=2\times 10^{5}$
subsequent MC iterations during which system properties were evaluated. The
$\hat{z}$ component of magnetization for a given sublattice $i$ (where the
number of sublattices is simply the number of magnetic ions in the magnetic
unit cell) is calculated as

$m_{z}^{i}=\frac{1}{N_{cells}}\sum_{j}^{N_{cells}}m_{z,j}^{i},$ (4)

where $N_{cells}$ is the number of magnetic cells in the MC simulation box
($42\times 42=1764$ and $32\times 19=722$ for the two surfaces studied).

##  Appendix C Multipolization tensors for different surfaces

Figure 5:  Depiction of unit cell bases used to calculated the multipolization tensor and associated surface magnetization for the various surfaces and methods. (a), (b), (c) and (d) correspond to first, second, third and four unit cells in Table 2 respectively. Green arrows depict the $3$ $\mathrm{\mu_{B}}$ $\mathrm{Cr^{3+}}$ magnetic moments, and the purple arrow indicates the formally $5$ $\mathrm{\mu_{B}}$ $\mathrm{Fe^{3+}}$ magnetic moment. (a) $(001)$ pristine $\mathrm{Cr_{2}O_{3}}$ surface, assuming the top-most $\mathrm{Cr}$ is magnetically ordered. (b) $(\bar{1}02)$ surface. (c) Realistic $(001)$ $\mathrm{Cr_{2}O_{3}}$ surface, discussed in the main text, where the topmost $\mathrm{Cr}$ layer is paramagnetic. Here, the multipolizaiton tensor is calculated using a basis which can be periodically tiled, and then no additional contribution is added from the paramagnetic layer, giving a total surface magnetization $\mu_{surf}^{total}$ which is the equal to $\mathcal{M}_{zz}$ from the bulk unit cell. (d) Unit cell corresponding to bulk periodic part of $(001)$ $\mathrm{Cr_{2}O_{3}}$, plus contribution from nonperiodic $\mathrm{Fe}$ and flipped $\mathrm{Cr}$ near the surface. The magnetic moments for the nonperiodic part divided by the cross-sectional area of the $(001)$ unit cell give $+9.56\mu_{B}/\mathrm{nm^{2}}$, which is added to the multipolization $\mathcal{M}_{zz}=-16.74\mu_{B}/\mathrm{nm^{2}}$ From the bulk unit cell to give a total estimated $\mu_{surf}^{total}$ of $+7.18\mu_{B}/\mathrm{nm^{2}}$ along $[001]$. Table 2: Positions, in fractional coordinates, and corresponding magnetic moments, all oriented along $[001]$, for magnetic ions in the four unit cells shown in Figure 5, used to calculate the multipolization tensors. Cartesian components, in $\mathrm{\AA}$, for the lattice vectors of the hexagonal $[001]$ oriented cells are $\mathbf{a}=[2.46,-4.26,0.0]$, $\mathbf{b}=[2.46,4.26,0.0]$ and $\mathbf{c}=[0,0,13.53]$. For the $[001]$ oriented surfaces, the moments only have a $\hat{z}$ component as $\hat{z}\parallel[001]$. The lattice vectors for the monoclinic $[\bar{1}02]$ oriented cell are $\mathbf{a}=[2.46,-4.26,0.0]$, $\mathbf{b}=[-2.46,1.42,-4.51]$ and $\mathbf{c}=[-4.92,2.84,4.51]$ with the standard $\hat{z}\parallel[001]$ Cartesian basis. They are $\mathbf{a}=[4.92,0.0,0.0]$, $\mathbf{b}=[0.0,5.33,0.0]$ and $\mathbf{c}=[0.0,-0.79,7.21]$ in the rotated Cartesian basis where $\hat{z}\parallel[\bar{1}02]$. The $\hat{x}$, $\hat{y}$, and $\hat{z}$ components of the $[001]$-oriented magnetic moments are given in this rotated basis. | ordered $(001)$ surface | ($(\bar{1}02)$ pristine surface | paramagnetic $(001)$ surface |  $(001)$ surface with $\mathrm{Fe}$ monolayer
---|---|---|---|---
$\mathrm{Cr}$ site | $r^{\mathbf{a}}$ | $r^{\mathbf{b}}$ | $r^{\mathbf{c}}$ | $\mu^{z}$ | $r^{\mathbf{a}}$ | $r^{\mathbf{b}}$ | $r^{\mathbf{c}}$ | $\mu^{x}$ | $\mu^{y}$ | $\mu^{z}$ | $r^{\mathbf{a}}$ | $r^{\mathbf{b}}$ | $r^{\mathbf{c}}$ | $\mu^{z}$ | $r^{\mathbf{a}}$ | $r^{\mathbf{b}}$ | $r^{\mathbf{c}}$ | $\mu^{z}$
1 | 0.33 | 0.67 | 0.014 | -3 | 0.694 | 0.389 | 0.056 | 0 | -2.54 | +1.60 | 0.67 | 0.33 | 0.00 | +3 | 0.33 | 0.67 | 0.00 | +3
2 | 0.00 | 0.00 | 0.153 | +3 | 0.194 | 0.889 | 0.056 | 0 | -2.54 | +1.60 | 0.33 | 0.67 | 0.028 | -3 | 0.67 | 0.33 | 0.167 | +3
3 | 0.67 | 0.33 | 0.181 | -3 | 0.194 | 0.50 | 0.25 | 0 | 2.54 | -1.60 | 0.00 | 0.00 | 0.167 | +3 | 0.33 | 0.67 | 0.195 | -3
4 | 0.33 | 0.67 | 0.319 | +3 | 0.694 | 0.00 | 0.25 | 0 | 2.54 | -1.60 | 0.67 | 0.33 | 0.195 | -3 | 0.00 | 0.00 | 0.33 | +3
5 | 0.00 | 0.00 | 0.347 | -3 | 0.694 | 0.389 | 0.56 | 0 | -2.54 | +1.60 | 0.33 | 0.67 | 0.333 | +3 | 0.67 | 0.33 | 0.361 | -3
6 | 0.67 | 0.33 | 0.486 | +3 | 0.194 | 0.889 | 0.56 | 0 | -2.54 | +1.60 | 0.00 | 0.00 | 0.361 | -3 | 0.33 | 0.67 | 0.50 | +3
7 | 0.33 | 0.67 | 0.514 | -3 | 0.194 | 0.50 | 0.75 | 0 | 2.54 | -1.60 | 0.67 | 0.33 | 0.50 | +3 | 0.00 | 0.00 | 0.528 | -3
8 | 0.0 | 0.0 | 0.653 | +3 | 0.694 | 0.00 | 0.75 | 0 | 2.54 | -1.60 | 0.33 | 0.67 | 0.528 | -3 | 0.67 | 0.33 | 0.667 | +3
9 | 0.67 | 0.33 | 0.681 | -3 | - | - | - | - | - | - | 0.00 | 0.00 | 0.667 | +3 | 0.33 | 0.67 | 0.695 | -3
10 | 0.33 | 0.67 | 0.819 | +3 | - | - | - | - | - | - | 0.67 | 0.33 | 0.695 | -3 | 0.00 | 0.00 | 0.833 | +3
11 | 0.00 | 0.00 | 0.847 | -3 | - | - | - | - | - | - | 0.33 | 0.67 | 0.833 | +3 | 0.67 | 0.33 | 0.861 | -3
12 | 0.67 | 0.33 | 0.986 | +3 | - | - | - | - | - | - | 0.00 | 0.00 | 0.861 | -3 | 0.00 | 0.00 | 1.028 | -3

In Table 2 We give the positions (in fractional coordinates) of the magnetic
ions, as well as the magnetic moments, in the unit cell bases which are used
to calculate the multipolization tensors corresponding to magnetization on the
$(001)$ and ($\bar{1}02)$ surfaces of $\mathrm{Cr_{2}O_{3}}$. Recall that with
a local moment approximation, the multipolization tensor $\mathcal{M}$ can be
calculated as:

$\mathcal{M}_{ij}^{lm}=\frac{1}{V}\sum_{\alpha}R_{i,\alpha}m_{j,\alpha},$ (5)

where $V$ is the unit cell volume, and the sum is over magnetic atoms in the
unit cell. Equation 5 is equivalent to Equation 1 in the main manuscript. By
inspection of the form of $\mathcal{M}$ it is clear that the
$ij^{\mathrm{th}}$ component of the tensor should ideally correspond to the
component of magnetization which is oriented along $\hat{i}$ on a surface
whose normal is parallel to $\hat{j}$. We remind the reader however that for a
given Miller plane with its nonpolar termination, only the out-of-plane
dimension of the corresponding bulk unit cell is unambiguously determined;
each in-plane lattice constant for the unit cell used to calculate
$\mathcal{M}$ can correspond to any arbitrary branch of multipolization
increment which is parallel to the surface normal. Therefore, to reliably
calculate the three cartesian components of magnetization on a given surface,
one should calculate $\mathcal{M}$ within a rotated basis where the $\hat{z}$
cartesian axis is parallel to the surface normal. This is already the case for
the $(001)$ surface with the conventional hexagonal unit cell. To obtain
$\mathcal{M}$ for the $(\bar{1}02)$ surface, we rewrite the lattice vectors in
a rotated cartesian coordinate system with $\hat{z}\parallel\hat{\mathbf{n}}$
where $\hat{\mathbf{n}}$ is the surface normal $[\bar{1},0,2]$, and $\hat{x}$
is parallel to the in-plane $\mathbf{a}$ lattice vector; The lattice vectors
in this rotated basis are also given in the table caption. The $\hat{x}$,
$\hat{y}$, and $\hat{z}$ components of the $[001]$-oriented $3\mu_{B}$
$\mathrm{Cr}$ magnetic moments are obtained simply by taking the dot product
of the $[001]$ direction with the rotated cartesian vectors. From this we
predict an in-plane (out-of-plane) $(\bar{1}02)$ surface magnetization of
$+7.55$ ($-4.75$) $\mu_{B}/\mathrm{nm^{2}}$ respectively as stated in the
manuscript. The in-plane magnetization is fully along the rotated cartesian
$\hat{y}$ direction corresponding to a nonzero $\mathcal{M}_{zy}$ component.
The corresponding bulk unit cell and basis for the $(\bar{1}02)$ surface is
also depicted in Figure 5(b).
We point out here that the surface magnetization is only rigorously tied to
the bulk multipolization tensor in the absence of any surface reconstruction,
spin disorder, spin flipping, or doping; if the surface is truly just an
abrupt termination of the bulk, the surface magnetization can then be
determined by calculating $\mathcal{M}$ for the bulk unit cell which
periodically tiles the semi-infinite surface of interest completely
analogously to the procedure for determining bound surface charge from
$\mathbf{P}_{\mathrm{bulk}}$. This is the case for the first two columns of
Table 2, representing respectively the $(001)$ surface assuming the surface
$\mathrm{Cr}$ retain the bulk AFM order (corresponding to the basis in Figure
5(a), identical to that in the left-hand side of Figure 1(b) in the main
text), and the $(\bar{1}02)$ surface. Recall from the main text that since
$\lambda_{surf}\sim\lambda_{bulk}$ for the $(\bar{1}02)$ Miller plane, the
outermost $\mathrm{Cr}$ have bulk-like magnetization and the mutlipolization
calculated from Equation 5 in this case corresponds rigorously to the bulk
$\mathcal{M}$ and should yield the true surface magnetization. However, for
the other two cases we discuss, i.e. realistic $(001)$ $\mathrm{Cr_{2}O_{3}}$
with a paramagnetic $\mathrm{Cr}$ surface layer, and $(001)$
$\mathrm{Cr_{2}O_{3}}$ with a single monolayer of $\mathrm{Fe}$ there is no
bulk unit cell which can be tiled semi-infinitely to define the surface.
One way to approximately calculate the surface magnetization in this case is
the method discussed in the main manuscript for $(001)$ $\mathrm{Cr_{2}O_{3}}$
with a paramagnetic surface layer. Here, a bulk multipolization tensor is
calculated based on the unit cell which can be tiled parallel to the surface
normal up to where the material maintains bulk character. For pristine $(001)$
$\mathrm{Cr_{2}O_{3}}$ this corresponds to unit cell shown in Figure 5(c),
identical to that on the right-hand side of Figure 1(b) in the main
manuscript. For $\mathrm{Fe}$-capped $\mathrm{Cr_{2}O_{3}}$ it corresponds to
the non-rectangular unit cell shown in Figure 5(d), which excludes the
$\mathrm{Fe}$ monolayer and the $\mathrm{Cr}$ which has flipped with respect
to the bulk AFM magnetic order. The corresponding positions, in direct
coordinates, for the magnetic $\mathrm{Cr}$ in these units cells are also
given in Table 2, yielding $\mathcal{M}_{zz}$ components of $-2.4$
$\mu_{B}/\mathrm{nm}^{2}$ and $-16.7$ $\mu_{B}/\mathrm{nm^{2}}$ respectively.
Then, to calculate the full theoretical magnetization for the actual surface,
one adds to the multipolization tensor-based value from this periodic unit
cell the remaining nonperiodic contribution. For Figure 5(c) the magnetization
contribution from the outermost paramagnetic layer is just zero, whereas for
Figure 5(d) this can be approximated by summing the magnetic moment vectors
for $\mathrm{Fe}$ and the flipped $\mathrm{Cr}$ ($5-3=+2$ $\mu_{B}$) and
dividing by the cross-sectional area of the $[001]$-oriented unit cell,
yielding $+9.56$ $\mu_{B}/\mathrm{nm^{2}}$. Adding these “nonperiodic”
contributions to the $M_{zz}$ components of the multipolization tensors as
calculated from equation 5 for the bulk unit cells gives $-2.4$
$\mu_{B}/\mathrm{nm^{2}}$ (as quoted in the manuscript) and
$-16.74+9.56=-7.18$ $\mu_{B}/\mathrm{nm^{2}}$ for $(001)$ pristine and
$\mathrm{Fe}$-capped surfaces respectively.

##  Appendix D Feasibility of synthesising $\mathrm{Cr_{2}O_{3}}$ capped with
an $\mathrm{Fe}$ monolayer

To assess the feasibility of terminating a $\mathrm{Cr_{2}O_{3}}$ slab with
$\mathrm{Fe}$, we have calculated the relative stability a relaxed
$\mathrm{Cr_{11}Fe_{1}O_{36}}$ slab structure with $\mathrm{Fe}$ replacing
$\mathrm{Cr}$ in the two positions directly below the topmost oxygen layer
($\mathrm{Cr}$ $11$ and $\mathrm{Cr}$ $10$ positions according to the labeling
in Figure 1(a) of the main text). We estimate the liklihood of substitution of
$\mathrm{Fe}$ on these sites rather than the terminating $12^{\mathrm{th}}$
position by first calculating the total energies within DFT+U for the spin-
polarized slabs with $\mathrm{Fe}$ in the $11^{\mathrm{th}}$ and
$\mathrm{10{th}}$ positions after fully relaxing the unit cell-thick slab. We
then calculate the $\mathrm{Fe}$-$\mathrm{Cr}$ couplings for the structure
with $\mathrm{Fe}$ in these intermediate positions via the usual total energy
method described earlier. We next calculate the Heisenberg contribution to the
total DFT+U total energy using these $\mathrm{Fe}$-$\mathrm{Cr}$ $J$ values
(as well as the $\mathrm{Cr}$-$\mathrm{Cr}$ $J$s for the $\mathrm{Cr}$ in the
unit cell) along with the relative directions of the magnetic moments in the
DFT+U calculation. Finally, we subtract off the magnetic Heisenberg
contribution from the DFT+U total energy; because the likelihood of site
substitution is primarily dependent on the atomic environment, we can get a
better idea of relative formation stability by neglecting magnetic
contributions to the energy. We find that the resulting energies for
$\mathrm{Fe}$ in the $11^{\mathrm{th}}$ and $10^{\mathrm{th}}$ positions are
$0.877$ $\mathrm{eV}$ and $0.896$ $\mathrm{eV}$ higher respectively than the
structure in Figure 4(c) of the main text with a terminating $\mathrm{Fe}$
layer in the $12^{\mathrm{th}}$ position. Thus at room temperature
$\mathrm{k_{B}T}\sim 0.025$ $\mathrm{eV}$, the probability of $\mathrm{Fe}$
substituting $\mathrm{Cr}$ at these sites is suppressed, and introducing
$\mathrm{Fe}$ into a vacuum chamber at the very end of growth should lead to a
reasonably uniform $\mathrm{Fe}$ monolayer on the top of
$\mathrm{Cr_{2}O_{3}}$. We recognize that synthesis of the final hypothetical
structure discussed in the main text, with a single $\mathrm{Fe}$ terminating
a $[001]$ oriented slab of $\mathrm{Cr_{2}O_{3}}$, is nontrivial, since the
oxygen termination of the $(001)$ $\mathrm{Cr_{2}O_{3}}$ structure which would
allow subsequent deposition of the $\mathrm{Fe}$ monolayer is polar.
Nevertheless, we believe it would be feasible, particularly given several
experimental studies indicated that the oxygen termination of chromia can be
stabilized by varying the oxygen partial pressure during synthesisLübbe and
Moritz (2009); Bikondoa  _et al._ (2010); Kaspar  _et al._ (2013).

## References

  * Astrov (1960) D. N. Astrov, J. Exptl. Theoret. Phys. (U.S.S.R.) 11, 984 (1960).
  * Belashchenko (2010) K. D. Belashchenko, Physical Review Letters 105 (2010), 10.1103/PhysRevLett.105.147204.
  * Ye (2022) S. Ye, Physica Status Solidi - Rapid Research Letters 16 (2022), 10.1002/pssr.202100396.
  * Hedrich  _et al._ (2021) N. Hedrich, K. Wagner, O. V. Pylypovskyi, B. J. Shields, T. Kosub, D. D. Sheka, D. Makarov, and P. Maletinsky, Nature Physics 17, 574 (2021).
  * Nogue´s and Schuller (1999) J. N. Nogue´s and I. K. Schuller, Journal of Magnetism and Magnetic Materials 192, 203 (1999).
  * Stamps (2000) R. L. Stamps, Journal of Physics D: Applied 33 (2000).
  * Echtenkamp (2021) W. Echtenkamp, “Voltage-controlled magnetization in chromia-based magnetic heterostructures heterostructures,”  (2021).
  * Schlitz  _et al._ (2018) R. Schlitz, T. Kosub, A. Thomas, S. Fabretti, K. Nielsch, D. Makarov, and S. T. Goennenwein, Applied Physics Letters 112 (2018), 10.1063/1.5019934.
  * Muduli  _et al._ (2021) P. Muduli, R. Schlitz, T. Kosub, R. Hübner, A. Erbe, D. Makarov, and S. T. Goennenwein, APL Materials 9 (2021), 10.1063/5.0037860.
  * Appel  _et al._ (2019) P. Appel, B. J. Shields, T. Kosub, N. Hedrich, R. Hübner, J. Faßbender, D. Makarov, and P. Maletinsky, Nano Letters 19, 1682 (2019).
  * Wörnle  _et al._ (2021) M. S. Wörnle, P. Welter, M. Giraldo, T. Lottermoser, M. Fiebig, P. Gambardella, and C. L. Degen, Physical Review B 103 (2021), 10.1103/PhysRevB.103.094426.
  * Spaldin (2021) N. A. Spaldin, Journal of Experimental and Theoretical Physics 132, 493 (2021).
  * Stengel (2011) M. Stengel, Physical Review B - Condensed Matter and Materials Physics 84 (2011), 10.1103/PhysRevB.84.205432.
  * Nakagawa  _et al._ (2006) N. Nakagawa, H. Y. Hwang, and D. A. Muller, Nature Materials 5, 204 (2006).
  * Vanderbilt and King-Smith (1993) D. Vanderbilt and R. D. King-Smith, Physical Review B 48, 4442 (1993).
  * King-Smith and Vanderbilt (1993) R. D. King-Smith and D. Vanderbilt, Physical Review B 47, 1651 (1993).
  * Spaldin  _et al._ (2013) N. A. Spaldin, M. Fechner, E. Bousquet, A. Balatsky, and L. Nordström, Physical Review B - Condensed Matter and Materials Physics 88 (2013), 10.1103/PhysRevB.88.094429.
  * He  _et al._ (2010) X. He, Y. Wang, N. Wu, A. N. Caruso, E. Vescovo, K. D. Belashchenko, P. A. Dowben, and C. Binek, Nature Materials 9, 579 (2010).
  * Fechner  _et al._ (2018) M. Fechner, A. Sukhov, L. Chotorlishvili, C. Kenel, J. Berakdar, and N. A. Spaldin, Physical Review Materials 2 (2018), 10.1103/PhysRevMaterials.2.064401.
  * Samuelsent  _et al._ (1970) E. J. Samuelsent, M. T. Hutchings, and G. Shirane, Physica 48 (1970).
  * Shi  _et al._ (2009) S. Shi, A. L. Wysocki, and K. D. Belashchenko, Physical Review B - Condensed Matter and Materials Physics 79 (2009), 10.1103/PhysRevB.79.104404.
  * Xiang  _et al._ (2011) H. J. Xiang, E. J. Kan, S. H. Wei, M. H. Whangbo, and X. G. Gong, Physical Review B - Condensed Matter and Materials Physics 84 (2011), 10.1103/PhysRevB.84.224429.
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Physical Review B 54 (1996).
  * Urru and Spaldin (2022) A. Urru and N. A. Spaldin, arXiv:2206.00522 (2022).
  * Mostovoy  _et al._ (2010) M. Mostovoy, A. Scaramucci, N. A. Spaldin, and K. T. Delaney, Physical Review Letters 105 (2010), 10.1103/PhysRevLett.105.087202.
  * Skubic  _et al._ (2008) B. Skubic, J. Hellsvik, L. Nordström, and O. Eriksson, Journal of Physics Condensed Matter 20 (2008), 10.1088/0953-8984/20/31/315203.
  * Wysocki  _et al._ (2012) A. L. Wysocki, S. Shi, and K. D. Belashchenko, Physical Review B - Condensed Matter and Materials Physics 86 (2012), 10.1103/PhysRevB.86.165443.
  * Kota  _et al._ (2013) Y. Kota, H. Imamura, and M. Sasaki, Applied Physics Express 6 (2013), 10.7567/APEX.6.113007.
  * Borisov  _et al._ (2005) P. Borisov, A. Hochstrat, X. Chen, W. Kleemann, and C. Binek, Physical Review Letters 94 (2005), 10.1103/PhysRevLett.94.117203.
  * Kota  _et al._ (2014) Y. Kota, H. Imamura, and M. Sasaki, IEEE Transactions on Magnetics 50 (2014), 10.1109/TMAG.2014.2324014.
  * Dedkov  _et al._ (2007) Y. S. Dedkov, C. Laubschat, S. Khmelevskyi, J. Redinger, P. Mohn, and M. Weinert, Physical Review Letters 99 (2007), 10.1103/PhysRevLett.99.047204.
  * Rosenberg and Franz (2012) G. Rosenberg and M. Franz, Physical Review B - Condensed Matter and Materials Physics 85 (2012), 10.1103/PhysRevB.85.195119.
  * Schulz  _et al._ (2019) S. Schulz, I. A. Nechaev, M. Güttler, G. Poelchen, A. Generalov, S. Danzenbächer, A. Chikina, S. Seiro, K. Kliemt, A. Y. Vyazovskaya, T. K. Kim, P. Dudin, E. V. Chulkov, C. Laubschat, E. E. Krasovskii, C. Geibel, C. Krellner, K. Kummer, and D. V. Vyalikh, npj Quantum Materials 4 (2019), 10.1038/s41535-019-0166-z.
  * Blochl (1994) P. E. Blochl, Physical Review B 50, 24 (1994).
  * Anisimov  _et al._ (1997) V. I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein, J. Phys.: Condens. Matter 9, 76603 (1997).
  * Dudarev  _et al._ (1998) S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Physical Review B 57 (1998).
  * Dudko  _et al._ (1971) K. L. Dudko, V. V. Eremenko, and L. M. Semenenko, physica status solidi (b) 43, 471 (1971).
  * Tobia  _et al._ (2010) D. Tobia, E. D. Biasi, M. Granada, H. E. Troiani, G. Zampieri, E. Winkler, and R. D. Zysler (2010).
  * Lübbe and Moritz (2009) M. Lübbe and W. Moritz, Journal of Physics Condensed Matter 21 (2009), 10.1088/0953-8984/21/13/134010.
  * Bikondoa  _et al._ (2010) O. Bikondoa, W. Moritz, X. Torrelles, H. J. Kim, G. Thornton, and R. Lindsay, Physical Review B - Condensed Matter and Materials Physics 81 (2010), 10.1103/PhysRevB.81.205439.
  * Kaspar  _et al._ (2013) T. C. Kaspar, S. E. Chamberlin, and S. A. Chambers, Surface Science 618, 159 (2013).