Surface Magnetization in Antiferromagnets: Classification, example materials, and relation to magnetoelectric responses

We use symmetry analysis and density functional theory to characterize antiferromagnetic (AFM) materials which have a finite equilibrium magnetization density on particular surface terminations. A nonzero magnetic dipole moment per unit area or"surface magnetization"can arise on particular surfaces of many AFMs due to the bulk magnetic symmetries. Such surface magnetization plays an essential role in numerous device applications, from random-access magnetoelectric (ME) memory to exchange bias. However, at this point a universal description of AFM surface magnetization is lacking. We first introduce a classification system based on whether the surface magnetization is sensitive or robust to roughness, and on whether the surface of interest is magnetically compensated or uncompensated in the bulk magnetic ground state. We show that uncompensated surface magnetization can be conveniently described in terms of ME multipoles at the local-moment, unit cell level, and demonstrate that the symmetry of the multivalued"multipolization lattice"distinguishes between roughness-robust and roughness-sensitive surface magnetization. We then demonstrate that magnetization on bulk-compensated surfaces arises due to ME multipoles (in addition to higher-order magnetic terms) at the atomic site level. These can further be understood in terms of bulk ME responses, arising from the effective electric field resulting from the surface termination. We also show with density functional calculations that nominally compensated surfaces in Cr2O3 and FeF2 develop a finite magnetization density at the surface, in agreement with our predictions based on both group theory and the linear and higher-order ME response tensors. Our analysis provides a comprehensive basis for understanding the surface magnetic properties in AFMs, and has important implications for phenomena such as exchange bias coupling.


I. INTRODUCTION
By definition, antiferromagnets (AFMs) have zero net magnetization, since the magnetic dipole moments of the oppositely oriented equivalent sublattices sum to zero over the bulk magnetic unit cell.However, theoretical arguments [1,2] and subsequent experimental measurements [3][4][5][6] have definitively demonstrated that AFMs possess a finite magnetic dipole moment per unit area on specifically oriented surface terminations due to the subset of bulk symmetries retained at the surface.We refer to this as "surface magnetization" throughout this article.Such surface magnetization enhances the utility of spintronicsbased devices that use bulk AFM domains as logical bits.Specifically, since the direction of surface magnetization is connected to the underlying bulk AFM domain, it is a directly detectable indicator of the domain state [7].Additionally, surface magnetization is believed to play a major role in exchange bias, whereby the magnetization direction of a ferromagnet (FM) is pinned by exchange coupling to an adjacent AFM [8].Exchange bias is exploited in a wide variety of applications, from magnetic random access memory cells to giant magnetoresistive read heads [9].
The most intuitive case of surface magnetization occurs when the vacuum termination, without further modification of the bulk AFM dipole order, leads to a layer of uncompensated magnetic dipoles at the surface.This is because the oppositely pointed dipole layer above it, which would compensate this layer in the bulk, gets cut off at a surface.We refer to this form of surface magnetization, depicted in Figs.1(a) and 1(b), as "uncompensated" in the manuscript.
But surface magnetization arising from uncompensated magnetic dipole moments is just one category.Indeed, Andreev pointed out nearly three decades ago that, for a general low-symmetry surface of any AFM (not just those with uncompensated moments), a finite surface magnetization is expected to develop [1].
This leads to the counterintuitive result that surface magnetization can also arise on crystallographic planes which have an AFM order in the bulk.While such order formally results in a vanishing magnetic dipole per unit area, symmetry reduction due to the surface termination allows for induced components of the magnetic moments, producing a net ferromagnetic (or ferrimagnetic) order and, thus, a nonzero surface magnetization.We refer to this as "induced" surface magnetization, depicted in Figs.1(c) and 1(d).
Furthermore, both uncompensated and induced surface magnetization can be distinguished based on whether the sign and magnitude of the surface magnetization are robust to macroscopic surface roughness, which we refer to as "roughness robust," or whether the surface magnetization switches sign between atomic steps, in which case we call it "roughness sensitive."As we discuss extensively later in this manuscript, the roughness-sensitive and roughnessrobust categories can be identified and distinguished for any surface orientation of interest based purely on the bulk and surface symmetries [10].Roughness-robust and roughness-sensitive uncompensated surface magnetization are depicted in Figs.1(a) and 1(b), and induced surface magnetization is shown in Figs.1(c) and 1(d) for roughnessrobust and roughness-sensitive cases, respectively.
The different categories of surface magnetization mentioned above have been discussed previously in various contexts [2,8,[11][12][13][14].Among them, a particularly relevant argument that lays the ground for the present work is the symmetry correspondence between roughness-robust surface magnetization and bulk magnetoelectric (ME) responses, that is, an induced bulk magnetization in response to an applied electric field in AFMs with broken time-reversal symmetry.This symmetry connection was first discussed by Krichevtsov et al. in the context of Kerr rotation measurements revealing nonreciprocal reflection of light [15].Later, Belashchenko motivated the correspondence from a fundamental standpoint by showing that the existence of a surface with normal n reduces the bulk point group symmetry in the same way as an electric field E [2].Thus, a vacuum-terminated surface in a ME AFM should develop the same magnetization components as those induced by an electric field applied to a bulk sample in the direction parallel to the surface normal [2,3].
To date, the ME community has focused primarily on the symmetry connection of roughness-robust surface magnetization to responses in linear ME AFMs.These are ME AFMs with broken time-reversal and inversion symmetries for which the strength of the induced bulk magnetization M is linearly proportional to the strength of the applied electric field E, M i ¼ P j α ij E j , where α is the linear ME response tensor [16].The prototypical linear ME AFM is Cr 2 O 3 (chromia), which has attracted interest for decades due to being the first linear ME AFM to be theoretically predicted and experimentally confirmed [16,17] as well as having an unusually high bulk AFM ordering temperature of approximately 300 K [18].The (001) surface of Cr 2 O 3 , i.e., the plane perpendicular to the staggered magnetization or Néel vector direction, has an experimentally confirmed out-ofplane, roughness-robust uncompensated surface magnetization.This has the same symmetry origin as the linear ME response in bulk Cr 2 O 3 along the [001] direction [4][5][6].Crucially, the direction of surface magnetization on (001) Cr 2 O 3 and other linear ME AFMs can be readily switched (along with switching of the underlying bulk AFM domain) by simultaneous application of electric and magnetic fields via so-called "magnetoelectric annealing" [3,9], thus rendering this type of surface magnetization relevant for controllable spintronics devices.
On the other hand, roughness-sensitive, uncompensated surface magnetization is frequently discussed qualitatively in the exchange-bias literature [8,11,19].In particular, the strength of exchange-bias-mediated magnetization pinning decreases with increasing roughness for most AFMs, and this is widely understood to be due to regions of the AFM surface with oppositely oriented signs of magnetization [8].However, this roughness-sensitive surface magnetization, which occurs in many AFMs with no symmetry-allowed ME response, is generally seen as completely disconnected from roughness-robust surface magnetization in ME AFMs, since an analogous symmetry-based description is not widely known.
Finally, only a few specific examples of induced surface magnetization, which we discuss in detail later in this manuscript, have been described in the literature [2,13,20].While the previously mentioned symmetry considerations [2] for roughness-robust, uncompensated surface magnetization can also identify the existence of induced, roughness-robust surface magnetization, this symmetry analysis does not distinguish between the two.Instead, whether surface magnetization is induced or uncompensated previously had to be inferred on a case-by-case basis from knowledge of the bulk Néel vector direction with respect to the surface.Moreover, while prior authors have shown that group theory can also explain roughness-robust induced surface magnetization in AFMs in which the linear ME effect is forbidden, specifically, in centrosymmetric FeF 2 [13], no explanation for the seeming lack of corresponding ME response has been proposed.Related to the prior paragraph, to our knowledge, the additional category of roughness-sensitive induced surface magnetization [Fig.1(d)] has not even been discussed.Lastly, while generic symmetry arguments and a few experiments point to the existence of induced surface magnetization, no ab initio calculations demonstrating its energetic stability exist to our knowledge.This article aims to fill these gaps in the understanding of surface magnetization.In doing so, we hope to stimulate a reexamination of the role of surface magnetization in AFM materials, in particular, for crystallographic planes for which it was previously neglected, as well as motivate new applications and detection methods.We achieve this by first demonstrating that both roughness-robust and roughnesssensitive surface magnetization can be described on an equal footing by a single, rigorous symmetry formalism for both uncompensated and induced cases.Second, we map the group-theory methodology to an alternative description based on bulk magnetic multipoles, which are established as ground-state symmetry indicators for bulk ME responses [21,22].This mapping establishes a bulk-boundary connection that allows us to predict the presence of surface magnetization based on knowledge of the bulk multipoles and demonstrates for the first time that surface magnetization always maps to ME responses and vice versa.Our final goal is to show with ab initio calculations for a few technologically relevant AFMs that induced surface magnetization which is symmetry allowed on nominally compensated surfaces is indeed energetically favorable and, thus, should occur in real materials.
In part I, we present a group-theory-based formalism to identify surfaces that have a nonvanishing equilibrium magnetization in bulk AFMs and to determine whether the magnetization is roughness robust or roughness sensitive.We build directly on the seminal result of Ref. [2], which showed that, for roughness-robust cases [23], one can identify whether the surface has an equilibrium magnetization by finding the symmetry operations in the bulk magnetic point group leaving the direction of the surface normal invariant and checking if this surface subgroup is compatible with ferromagnetism [2].We extend this result and develop a concrete step-by-step procedure which considers the bulk magnetic space group, including translations, instead of the magnetic point group, thereby allowing us to identify and distinguish roughnessrobust and roughness-sensitive surface magnetization.
In parts II and III, we connect the categories of surface magnetization to the previously mentioned bulk magnetic multipoles, which capture higher-order asymmetries beyond magnetic dipole order in a generic, spatially inhomogenous bulk magnetization density μðrÞ [24][25][26].In part II, we begin with cases for which the bulk magnetic multipoles corresponding to the surface magnetization are the ME multipoles, M ij ¼ R r i μ j ðrÞd 3 r, with i and j referring to Cartesian components (the ME multipoles are also referred to as "magnetic quadrupoles" in the literature [27]).These form the next-order term in the multipole expansion after the magnetic dipoles and have the same symmetry as the linear ME effect.Then, we show that the dominant contribution to uncompensated surface magnetization can always be captured by a simplified expression, the so-called "local-moment" contribution, that treats the ME multipole as a sum over local dipole moments centered at the atomic positions in the unit cell, without taking into account the full inhomogeneity of the magnetization density μðrÞ [28].We find that the connection of uncompensated surface magnetization to the bulk local-moment ME multipole per unit volume (the ME multipolization) is analogous to the link between surface charge and bulk electrostatic polarization [29].In particular, in the same way that the polarization is well established to be multivalued within the modern theory of polarization [30], the local-moment part of the ME multipolization is multivalued due to the periodicity of the lattice [24,28].This multivaluedness, in turn, allows us to distinguish roughnessrobust from roughness-sensitive uncompensated surface magnetization, with roughness-sensitive surface magnetization always corresponding to a centrosymmetric array of allowed ME multipolization values.We illustrate both cases, choosing (001) Cr 2 O 3 for the roughness-robust, uncompensated example, motivated by its previously mentioned relevance in spintronic applications, and compare it with roughness-sensitive (001) Fe 2 O 3 .Fe 2 O 3 is isostructural to Cr 2 O 3 but due to its inversion-preserving magnetic order has no macroscopic linear ME response.We also work through the example of roughness-sensitive uncompensated surface magnetization in (111) NiO, which has promise in a variety of energy-based applications [31].We provide a flowchart summarizing the step-by-step procedure to identify and categorize uncompensated surface magnetization at the end of part II.
In part III, we move to induced surface magnetization and show that, while it is missed by the local dipole moment simplification used in part II, it is contained in the atomic-site multipoles that capture asymmetries of the magnetization density around the individual magnetic atoms [26,32,33].We extend the discussion of multipoles beyond the linear ME multipoles to second-and third-order multipoles, linked to quadratic and cubic ME responses.In doing so, we reveal the crucial point that, while uncompensated surface magnetization can always be described by the local-moment ME multipoles, the atomic-site multipoles which indicate induced surface magnetization do not necessarily have to correspond to a linear ME effect.Instead, some cases of surface magnetization are described by higher-order atomic-site multipoles which correspond to local ME responses that may be quadratic, cubic, etc., in the strength of the electric field [22,32].This is because the group-theory prescriptions to identify surface magnetization are equivalent to the formalism to identify the ME responses in the bulk at all orders in the electric field's strength E. This point has to date been overlooked.
In Sec.IV B, we analyze roughness-robust, induced surface magnetization and show that it corresponds to atomic-site multipoles which are ferroically ordered throughout the bulk unit cell, implying corresponding net ME responses which are either linear, quadratic, or cubic in the electric field strength, depending on the symmetry of the multipole.We take the nominally compensated Cr 2 O 3 ð 120Þ and (100) and FeF 2 (110) surfaces as examples for which we confirm our multipole-based predictions of induced surface magnetization with ab initio calculations.Our choice of (110) FeF 2 as an additional example material is motivated by the experimentally confirmed, large exchange bias on this surface despite its nominally compensated surface magnetization [8].In general, the overlooked ubiquity of induced surface magnetization implies that exchange bias for AFM interfaces with compensated magnetic dipoles, which has confounded the community for decades, must be reexamined.
In Sec.IV C, we discuss how roughness-sensitive, induced surface magnetization arises due to each surface plane having ferroically ordered atomic-site multipoles which change sign on adjacent planes.Finally, in Sec.V, we examine the "null" case, i.e., compensated surfaces which do not have any of the four forms of surface magnetization in Fig. 1.We take as an example the ð1 10Þ surface of NiO and show with symmetry arguments and density-functional theory (DFT) calculations that this surface has zero symmetry-allowed uncompensated or induced magnetization due to the in-plane antiferroic ordering of both its dipoles and all of its atomic-site multipoles.We also provide a flowchart at the end of part III depicting the procedure to identify and classify induced magnetization on surfaces with no uncompensated magnetization.Section VI contains our conclusions.
In addition to providing a unifying symmetry formalism that synthesizes multiple distinct literature reports on AFM surface magnetization [1,2,13,34,35], our work introduces new perspectives for understanding, as well as exploiting, AFM surface magnetization.First, it demonstrates that all categories of surface magnetization in AFMs can be conveniently described in terms of bulk magnetic multipoles, implying that bulk multipoles may serve as order parameters for surface magnetization, and, conversely, surface magnetization may indicate the existence of hidden bulk multipoles.Moreover, the universal connection of surface magnetization to magnetic multipoles demonstrates that all categories of surface magnetization correspond to ME responses, which in the roughness-sensitive cases cancel at the bulk level and in the roughness-robust cases correspond to a finite bulk ME effect.This suggests that AFM surface magnetization may be effectively controlled via ME annealing in many more materials than was previously thought.Overall, our findings show that the phase space of bulk AFMs and crystallographic planes with surface magnetization that can be exploited in applications is vastly richer than the cases currently discussed in the literature.

II. PART I: GROUP-THEORY DESCRIPTION AND CLASSIFICATION OF SURFACE MAGNETIZATION IN ANTIFERROMAGNETS
To start, we introduce the group-theoretical procedure we use to determine first whether a certain Miller plane for a particular AFM will allow for a finite equilibrium surface magnetization and second whether the surface magnetization is roughness robust or roughness sensitive.We emphasize that the group-theory procedure to identify roughness-robust surface magnetization has already been developed in Refs.[2,13].Here, we extend the procedure to include, and distinguish, the roughness-sensitive categories of surface magnetism [1].

A. Surface magnetization in the atomically flat limit
Let us consider a generic bulk AFM, periodic in all spatial directions, with magnetic space group (MSG) G ¼ fðRjtÞg.Here, fRg is the set of point group operations in G (some of which may contain time-reversal Θ), and ftg the set of fractional and Bravais lattice translations which, when combined with fRg, comprise the complete magnetic space group.We note that the prescriptions discussed below apply in the same way to collinear as well as noncollinear AFMs; in the examples, we consider collinear AFMs for simplicity.Now, we take a particular Miller plane ðhklÞ and create a surface by terminating the solid and interfacing it with vacuum.Assuming a macroscopically flat termination with no steps (we later relax this requirement), the unit vector surface normal n⊥ðhklÞ, which is perpendicular to the vacuum-terminated surface, is all that is needed to characterize the subgroup of the bulk MSG that defines the two-dimensional surface in question.This "surface" MSG, which we denote G n, is the set of space group operations fðR njt ⊥ nÞg in G that leave the polar vector n invariant [36], modulo translations parallel to the surface or, equivalently, perpendicular to n: where ðR i jt i Þ is the ith element of the bulk MSG G and ∧ is the logical and.Note that the surface MSG must have broken inversion symmetry I, since I always reverses the direction of n.Other operations in G which are preserved or broken in G n depend on the relative orientations of n and the principle axes of the bulk MSG operations.
Once we have identified the surface MSG G n for the Miller plane of interest, we next consider the corresponding magnetic point group (MPG) g n ¼ fR ng which contains the point group elements of G n without the accompanying translations.We then check whether g n, which is a subgroup of the bulk MPG, allows for ferromagnetism.31 out of the 58 MPGs in solids are FM compatible, meaning that at least one component of magnetization is left invariant under all operations of the point group [37].In three dimensions, a FM-compatible MPG implies that a finite magnetic dipole moment per unit volume is allowed.In the present formalism, g n describes the two-dimensional surface perpendicular to n: in this case, a FM-compatible surface MPG allows for a nonzero magnetic dipole moment per unit area, i.e., a surface magnetization.Note that the bulk three-dimensional AFM still has zero net magnetization per unit volume.Conversely, if g n is not FM compatible, an equilibrium surface magnetization on ðhklÞ⊥ n is symmetry forbidden.This simple procedure to identify atomically smooth surfaces of AFMs which are expected to have nonzero surface magnetization is illustrated with cartoons in Figs.2(a)-2(c).
As an aside, we mention that the procedure to identify g n starting from the bulk MSG relies on the absence of symmetry-breaking surface reconstructions or spin reorientations at the surface, which may occur for an experimental sample.In the case of surface reconstruction, it is still possible to identify the lower-symmetry, reconstructed surface MPG by finding the subset of operations in g n which leave the direction of n invariant and also leave invariant the positions of the atoms at the reconstructed surface.Since this MPG is necessarily a subgroup of g n, we can infer that if g n is FM compatible, then the reconstructed MPG must be as well.However, the specific allowed components of magnetization cannot be directly identified from the bulk surface MSG via the requirements in Eqs.(1).Similarly, surface spin reorientations or changes in the direction of the magnetic easy axis can change the surface MPG compared to the one obtained with the bulk magnetic ordering.In this case, it is possible to correctly identify the surface MPG by applying the procedure described earlier for pristine, unreconstructed surfaces, starting from the MSG of a bulk system with magnetic ordering reoriented according to the surface of interest.
In what follows, we explain how to determine whether the surface magnetization identified in Fig. 2 is robust in the presence of atomic steps or whether it averages to zero if roughness is introduced.

B. "Roughness-robust" versus "roughness-sensitive"
surface magnetization In the previous section, by restricting the surface MSG to operations which involve only translations parallel to the surface ⊥ n, we show how to identify AFM surfaces which have a macroscopic magnetization in the limit of an ideally smooth plane with no roughness or steps.Next, we relax this restriction and consider symmetry operations which connect flat surface regions separated by translations with a component perpendicular to the surface.This allows us to distinguish between FM-compatible surfaces whose macroscopic magnetization is robust to roughness (corresponding Group-theory procedure to identify whether a given Miller plane ðhklÞ⊥ n has a symmetry-allowed surface magnetization, assuming that the surface has no atomic steps.(a) The bulk magnetic space group (MSG) G is identified for the bulk AFM which is periodic in all directions.(b) The surface MSG G n is identified.(c) If the magnetic point group (MPG) g n corresponding to G n is compatible with ferromagnetism, then ðhklÞ⊥ n has an equilibrium surface magnetization m surf ðg nÞ.
to the subset which Refs. [2,13] identify) and those whose magnetization is roughness sensitive.Assume that for some bulk AFM we have already identified a particular Miller plane characterized by surface normal n which has equilibrium surface magnetization according to the procedure in Fig. 2. Now, let us go back and again consider the higher-symmetry bulk MSG G.We want to check whether there exists in G at least one operation and Here, m FM ðg nÞ is the magnetization which remains invariant under the operations of the surface MPG g n.If Eqs. ( 2)-(4) all hold, then there exist atomic steps on surface ðhklÞ⊥ n, connected by fractional translation t i , which are energetically and symmetrically equivalent but which have opposite signs of surface magnetization [1,2].Thus, for a macroscopic region with surface roughness in thermodynamic equilibrium, steps connected by t i with opposite directions and equal magnitudes of magnetization occur with equal probability such that the total magnetization averaged over the entire surface area is zero [38].
Conversely, if all bulk MSG operations which involve t= ⊥ n leave m g n FM unchanged, then all atomic steps in the presence of roughness have the same direction of surface magnetization, and, hence, the surface magnetization is roughness robust.This is illustrated in Fig. 3.
The group-theory formalism which we have described thus far is general, identifying all AFM Miller planes having a symmetry-allowed finite equilibrium surface magnetization which either may exist despite roughness or may require an atomically flat surface to be realized.Figure 4 gives a table categorizing these four types of surface magnetization and indicating which are allowed, given the absence or presence of inversion and global timereversal symmetries (as well as combined time-reversal plus translation symmetries) in the bulk MSG.The importance of the inversion symmetry requirements, which we have not discussed yet, becomes clear when we introduce the multipole-based description of surface magnetization in parts II and III.We can see that the roughness-robust, uncompensated surface magnetization has the strictest symmetry requirements for the bulk AFM, namely, having a bulk MSG with broken time-reversal and inversion symmetries.Note that these are the same symmetries which must be broken for a linear ME response.On the other hand, the roughness-sensitive categories, both induced and uncompensated, can occur for bulk AFMs with all possible combinations of inversion and time-reversal symmetries.
In the following parts of the manuscript, we introduce a complementary description of AFM surface magnetization, combined with concrete material examples, based on the multipoles of the bulk magnetization density, which help us to classify and understand each of the four categories.We see that this multipole-based formalism is intimately connected to the group-theory description we have reviewed and modified above.At the same time, it provides new insights via explicitly connecting surface magnetization to bulk properties and corresponding bulk ME responses.

III. PART II: UNCOMPENSATED SURFACE MAGNETIZATION IN TERMS OF LOCAL-MOMENT MAGNETOELECTRIC MULTIPOLIZATION
In part II of this manuscript, we develop a semiquantitative description of uncompensated AFM surface magnetization (both roughness-robust and roughness-sensitive) in terms of a "local-moment" ME multipole tensor which involves summing over the product of positions and dipole moments of magnetic atoms in the bulk unit cell tiling the slab containing the surface in question.This "multipolization" formalism was already proposed in Ref. [28] to predict the value of uncompensated, roughness-robust surface magnetization in linear MEs.Here, we show that this local-moment multipolization can also describe and distinguish roughness-sensitive uncompensated surface magnetization, even when the bulk MSG forbids a net ME response.We first review the description of bulk magnetization density in terms of a multipole expansion and then investigate the implications of the form of the local-moment ME multipolization on uncompensated surface magnetization.

A. Magnetoelectric multipoles as bulk indicators of surface magnetization
While in many cases one can approximate the magnetic order of a solid as a series of point magnetic dipole moments localized on the magnetic ions, the full magnetization density is more complex.It is convenient to describe the asymmetry beyond the usual magnetic dipole of a generic, inhomogenous magnetization density in terms of a multipole expansion.Usually, one constructs the form of the magnetic multipoles at each order via an expansion of the interaction energy of a magnetization density μðrÞ interacting with a spatially varying magnetic field HðrÞ [28,39]: where ði; jÞ ¼ 1, 2, 3 refer to Cartesian components and summation over repeated indices is implied.The first term of the expansion in Eq. ( 5) is the usual Zeeman term, containing the magnetic dipole, whereas the coefficient of the second term, is known as the ME multipole tensor and describes the first order of asymmetry in magnetization density [24][25][26][27].We briefly mention here that the rank-3 tensor which gives the coefficient of the third term in Eq. ( 5), is called the magnetic octupole [22,39].Note that Eq. ( 7) is even under inversion symmetry I, and in part III we show that ferroically ordered octupoles explain the existence of roughness-robust, induced surface magnetization in centrosymmetric AFMs.For uncompensated surface magnetization, however, the lower-order ME multipoles are sufficient to describe both roughness-sensitive and roughness-robust cases, so we concentrate exclusively on M ij in this section.
It is evident that the ME multipole tensor Eq. ( 6) simultaneously breaks inversion symmetry and time reversal due to its dependence on both position and magnetization.In fact, one can show rigorously that nonzero values of the M ij tensor imply a corresponding nonzero linear ME response α ij in the presence of an applied electric field [25], where as a reminder the linear ME response tensor α is defined by with E the applied electric field and M the induced magnetization in the bulk.Thus, the linear ME effect, which is a thermodynamic response, can be conveniently understood in terms of the ME multipole tensor, which characterizes the second-order asymmetry in bulk magnetization density in the absence of applied fields.The symmetry connection between surface magnetization and bulk ME effects mentioned in the introduction already provides a hint that the ME multipole tensor may also be relevant in describing surface magnetization.Crucially for our discussion of surface magnetization categories, M can be decomposed into a sum of "local-moment" and "atomic-site" terms in the following way [21,40]: Here, R i;α is the ith Cartesian component of the position of atom α in the magnetic unit cell, m j;α is the jth Cartesian component of the magnetic dipole moment of atom α, the sum is over magnetic atoms in the unit cell, and Ω α is a spherical region around atom α.Note that Eq. ( 9) neglects contributions to M ij due to magnetization density in interstitial regions, which we assume to be negligible.The first term in the expression, the "atomic-site" (AS) term, originates from asymmetry in the magnetization density around an atomic site, and the second term, called the "local-moment" (LM) term, describes contributions to M due to the asymmetric arrangement of local dipole moments across all magnetic sites at the unit-cell level.The two contributions to the ME multipole tensor are depicted with a cartoon in Fig. 5.When considering M LM ij in an extended solid, we normalize the tensor by the volume V of the magnetic unit cell: We refer to the volume-normalized, local-moment contribution to the ME multipole tensor as the "ME multipolization."In part III, we show in detail that the atomic-site contribution to the ME multipolization, as well as higher-order atomic-site terms in the multipole expansion [22,39] of Eq. ( 5), are bulk indicators for the existence of induced surface magnetization.In this part, however, we discuss only the local-moment contribution given by Eq. (10), which, as we see, is sufficient to describe all categories of uncompensated surface magnetization in AFMs.
It is also important to note that, in a bulk periodic solid, there is freedom in selecting the basis for our magnetic unit cell; specifically, we can translate any magnetic moment m α by a Bravais lattice vector R without changing our system.Therefore, MLM ij is defined only modulo an "ME multipolization increment" corresponding to such a translation [24][25][26]28].In fact, since each of the l magnetic atoms in the unit cell can be translated by any Bravais lattice vector, there can be up to 3l linearly independent multipolization increments [24].Once we have made a concrete basis choice, however, from Eq. (10), it is clear that MLM ij calculated for a unit cell with AFM-coupled layers of inplane FM magnetic dipole moments stacked along the ith Cartesian direction must be nonzero.We now motivate the implications of this crucial fact on uncompensated surface magnetization by a brief comparison with the case of charge on the surfaces of materials with nonzero electric polarization.
Recall that there is a correspondence between the bound charge on a given surface σ surf and the component of the bulk electric polarization vector, P bulk perpendicular to the surface, that is, P bulk • n ¼ σ surf [30].However, within the modern theory of polarization, the periodicity of a bulk crystal implies that P bulk , analogously to the multivalued MLM ij , is defined only modulo a "polarization quantum" corresponding to translating one electron by a lattice vector [41].As a consequence, rather than a single value of P bulk , a given bulk crystal has a polarization "lattice" with values separated by the polarization quantum P q ¼ eR=V, with R a Bravais lattice vector.Therefore, at first glance, the correspondence to surface charge for a given termination P bulk • n ¼ σ surf seems ambiguous.This apparent issue is solved by the fact that 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 [29].Thus, the component of P bulk along the surface normal, P bulk • n, calculated for this particular unit-cell choice is single valued, and the surface charge is determined unambiguously.For a given Miller plane ðhklÞ, the atomic

ME multipole
Atomic-site Local-moment FIG. 5. Local-moment and atomic-site contributions to magnetic multipoles.Each order of magnetic multipoles beyond the magnetic dipoles can be decomposed into a sum of atomic-site terms, which capture asymmetries of μðrÞ around the atomic nuclei (blue circles), and a local-moment contribution due to the asymmetric spatial arrangement of magnetic dipoles (purple arrows) centered on the ions in the unit cell.We depict only the second-order ME multipoles here for convenience.
termination which is electrostatically stable is defined by the bulk unit cell for which P bulk • n ¼ σ surf ¼ 0 [42].
While P bulk has units of electric charge per unit area, the components of the volume-normalized M has units of magnetic dipole moment per unit area.Thus, just as the first-order term (polarization) in a multipole expansion of the interaction energy of a bulk charge density in an electric field gives rise to a zeroth-order term (charge) associated with its surface, in Ref. [28], we argue that the second-order term (the ME multipolization) in the multipole expansion in Eq. ( 5) should lead to a first-order term (magnetic dipole per unit area) associated with a surface.Specifically, the bulk multipolization component Mij should correspond to the magnetic dipole density per unit area along Cartesian direction | with surface unit normal along Cartesian direction { [43].This surface magnetization can, in general, have contributions from both the atomic-site and localmoment terms in the Mij tensor.As we pointed out before, the local-moment contribution to Mij , like P bulk , is multivalued for a fully periodic solid.But, again, the choice of a specific Miller plane and atomic termination fixes the basis of the bulk unit cell used to compute Eq. ( 10) along the direction of the surface normal.Thus, the three components j ¼ 1, 2, 3 of surface magnetization associated with Mij coming from local-moment and atomic-site contributions on surface ðhklÞ⊥ î are all single valued.
However, for a surface with uncompensated magnetic dipoles ferromagnetically ordered in the surface plane, the major contribution to the surface magnetization is quantitatively captured by the local-moment term MLM ij alone.We can see this by noting first that MLM ij is necessarily nonzero in such a situation, as we mentioned previously.Second, the local-moment term Eq. ( 10) weights the contribution of each magnetic moment in the unit cell by its distance from the surface, thus capturing the quantitative uncompensated surface magnetization.
We illustrate the correspondence between uncompensated surface magnetization and local-moment ME multipolization concretely in Fig. 6(a) for (001) Cr 2 O 3 .The top left in the figure shows the 30-atom bulk unit cell in the hexagonal setting, with an "up-down-up-down" AFM order of the Cr moments along [001].This surface defines the electrostatically stable surface termination [top right in Fig. 6(a)] of (001) chromia, with a single Cr ion above the final oxygen layer.Using the Cr positions in the bulk unit cell in a Cartesian basis where zk½001 and assuming the idealized 0 K limit where the Cr moments have their formal 3μ B value and are polarized fully along [001], Eq. ( 10) yields a z=½001 oriented magnetization of 11.9μ B =nm 2 on the (001) surface corresponding to MLM zz [28,44].On the other hand, we can instead choose the bulk unit cell on the left in Fig. 6(b), which can be obtained from the unit cell in Fig. 6(a) by translating the topmost Cr ion downward along ½00 1 by one lattice vector.The bulk unit cell in Fig. 6(b) defines the polar, oxygen-terminated (001) surface on the right in Fig. 6.Recalculating MLM zz using the Cr ion positions of this new unit cell yields a magnetization of −2.4μ B =nm 2 .The difference between these two values, 11.9 − ð−2.4Þ ¼ 14.3μ B =nm 2 , is precisely equal to the multipolization increment for the ðz; zÞ component of MLM , where jcj is the length of the hexagonal lattice vector along [001].We must caution that the rigorous quantitative correspondence of MLM ij to the surface magnetization values relies first on fully ordered surface magnetic moments, which, as we discuss in detail in Ref. [44], is often an inaccurate assumption at elevated temperatures, even if the bulk sublattices are fully ordered.Second, the quantitative accuracy of MLM ij assumes that surface reconstructions are absent.And, finally, the quantitative rigor of local-moment multipolization relies on the absence of changes in anisotropy or spin reorientation at the surface.We refer the reader to the supplemental material in Ref. [44] for a discussion of approximate methods to predict quantitative surface magnetization by combining the contributions from MLM ij with nonperiodic contributions due to temperature effects, reconstructions, or spin reorientations.However, the discussion is beyond the scope of this work, where we are concerned primarily with qualitative identification of categories of surface magnetization based on bulk order parameters.Now, we point out a salient feature of the so-called "multipolization increments."Like the lattice of electrical polarization values, the complete set of local-moment ME multipolization values corresponding to a surface with normal n must be invariant under all symmetry transformations in the MSG of the bulk AFM [24,45].Therefore, if the bulk AFM has inversion symmetry I, the array of allowed values for a given component MLM ij must take the following centrosymmetric form: where with a ki ¼ a k • î the projection of the kth Bravais lattice basis vector a k onto Cartesian direction î and Δ MLM kj ¼ m j ja k j=V is the generalized form of the ME multipolization increment Eq. (11).ξ k in the first term of Eq. ( 12) can take the values 0 or 1.In case ξ k ¼ 0 for all k, we refer to Eq. ( 12) as a "zero-containing" multipolization array, or a "half-increment-containing" multipolization array otherwise.In fact, because the local-moment ME multipolization in Eq. ( 10) involves a summation over only the magnetic atoms in the unit cell, multipolization increments for a particular material, and the corresponding allowed surface magnetization values, must take the form of Eq. ( 12) as long as the lattice occupied by the magnetic atoms is centrosymmetric, even if the bulk AFM including all magnetic and nonmagnetic atoms breaks I.We note that in most cases (including all the examples considered in this manuscript) one can choose the unit cell semi-infinitely tiling the surface of interest such that one lattice vector a k is parallel to the surface normal and the other two lattice vectors lie in the surface plane.Then, if we take a Cartesian basis with Cartesian direction îka k k n, it follows that Δ MLM ij;k ¼ Δ MLM kj , and the sum over k in Eq. ( 12) disappears for the multipolization array corresponding to magnetization on the surface perpendicular to î [46].
The centrosymmetric form of Eq. ( 12) implies that, for an AFM having such a multipolization array, a Miller plane with surface normal nk î that has a specific atomic termination corresponding to a ĵ-oriented surface magnetization μ j ð nk îÞ ¼ MLM ij also has an atomic termination having surface magnetization with opposite direction, −μ j ð nk îÞ ¼ − MLM ij .Moreover, for AFMs with bulk I symmetry, these two terminations must be symmetrically equivalent and, thus, energetically degenerate.We can conclude then that nominally uncompensated surface magnetization in bulk AFMs that are centrosymmetric must be roughness sensitive.We emphasize again, however, that roughness-robust induced surface magnetization is still possible in centrosymmetric AFMs such as FeF 2 ; we return to this case in Sec.IV B. This multipolization-based argument is, in fact, just a complementary way of viewing Eqs. ( 2)-( 4) in Sec.II B, which characterize roughnesssensitive surface magnetization for uncompensated surfaces within our group-theory formalism.
Although uncompensated surface magnetization is always roughness sensitive for I-symmetric AFMs, this category can also occur in AFMs with broken I for specific cases.One such example can occur when an I-broken AFM has a centrosymmetric magnetic lattice, as pointed out previously.A second possibility is that, even if the magnetic lattice breaks I symmetry, additional bulk MSG symmetries constrain the ME multipolization array to the centrosymmetric forms Eq. ( 12) along specific directions.This can occur, for example, if there exists a mirror plane containing the surface normal, and the magnetization is oriented parallel to the mirror plane.In both of these cases, provided there exists a bulk symmetry in the bulk MSG that connects the terminations with opposite signs of surface magnetization, the surface magnetization will be roughness sensitive [47].
Finally, AFMs with broken I in the bulk and a noncentrosymmetric magnetic lattice can have I-broken multipolization arrays for some directions, where MLM;0 ij is the "spontaneous" ME multipolization that occurs for AFMs with broken I and Θ symmetries [24,28] (analogously to the spontaneous electric polarization in ferroelectrics) and is neither zero nor half a multipolization increment.For the corresponding surfaces, a termination with surface magnetization μ j ð nk îÞ does not have a symmetry-connected, energetically degenerate termination with surface magnetization −μ j ð nk îÞ.Therefore, unlike inversion-symmetric AFMs, uncompensated surface magnetization in inversion-asymmetric AFMs can be roughness robust.
To conclude Sec.III A of this paper, in Fig. 7 we provide a flowchart with steps to categorize the primary form of surface magnetization for a bulk AFM and a given surface of interest (recall that, in general, induced surface magnetization can coexist with uncompensated surface magnetization, though in these cases the latter dominates).First, one finds the bulk unit cell which tiles semi-infinitely the desired surface.Second, one computes the local-moment ME multipolization MLM ij for a specific atomic termination, as well as the multipolization increment Δ MLM ij corresponding to the surface.If the multipolization array is noncentrosymmetric, the surface magnetization must be uncompensated and roughness robust.If the multipolization array is centrosymmetric and the bulk MSG contains I, the surface magnetization must be roughness sensitive.If the multipolization array is centrosymmetric, but the bulk MSG breaks I, one has to resort to the group-theory procedure in part I to distinguish between roughness-sensitive and roughness-robust uncompensated cases.
Finally, if MLM ij and Δ MLM ij are identically zero, the only possible categories of surface magnetization are induced.In order to categorize and understand these induced varieties of surface magnetization (as well as the "null" case where no surface magnetization is allowed), we need to analyze the atomic-site multipoles as opposed to the local-moment contribution.We describe this in detail in part III.FIG. 7. Flowchart describing the steps to identify, given a bulk AFM and the surface plane ðhklÞ⊥ n of interest, whether the surface has nonzero surface magnetization and, if so, the category (uncompensated, roughness-robust; uncompensated, roughness-sensitive, or induced).

B. Examples of uncompensated surface magnetization
1. Uncompensated, roughness-robust: (001) Let us first recall the minimal symmetry requirements which we give in Fig. 4 for the bulk AFM MSG in order to have at least one Miller plane with this type of surface magnetization.First, as we already discussed extensively in Sec.III A, uncompensated, roughness-robust surface magnetization can occur only in AFMs with broken inversion in the bulk MSG.Next, we consider the effect of time reversal Θ.For any surface normal n and any magnetization pseudovector m, ð n; mÞ → ð n; −mÞ under the action of Θ.Therefore, either any possible Miller plane will be incompatible with surface magnetization, or the surface magnetization has to be roughness sensitive based on Eqs. ( 2)-(4).It is, thus, impossible to have a surface with roughness-robust surface magnetization if Θ or ðΘjtÞ is a bulk symmetry.
The combined requirement of broken I and Θ implies that this type of surface magnetization can occur only in AFMs with a bulk MSG that allows for the linear ME effect [48].
We now work through the full MSG symmetry analysis for the prototypical case of (001) Cr 2 O 3 which we already briefly discussed in connection to our discussion of localmoment ME multipolization in Sec.III A. In Fig. 8, we show the structure of Cr 2 O 3 , both the primitive rhombohedral cell and the top view of the conventional hexagonal cell.The bulk MSG of linear ME Cr 2 O 3 is R 30 c 0 [167.106][49].Although the crystal lattice of Cr 2 O 3 is centrosymmetric, the AFM ordering breaks inversion symmetry in the MSG.Therefore, R 30 c 0 has neither space inversion nor time reversal as a symmetry, so we know immediately that it can host uncompensated, roughness-insensitive surface magnetization for at least some Miller planes.It contains clockwise and counterclockwise threefold rotations and rotoinversions about the [001] axis and twofold rotations with corresponding mirror planes perpendicular to the rotation axes along [100], [110], and [010].Because R 30 c 0 is a "type-III" MSG [49,50], half of the 36 space group operations contain time reversal.Indeed, although I and Θ are individually broken, the product IΘ is a symmetry of the bulk MSG.Now consider the unit vector nk½001 perpendicular to the (001) surface.The operations in R 30 c 0 that either leave n invariant modulo translations ⊥½001 are ðEj0Þ, ð3 þ ½001 j0Þ, and ð3 − ½001 j0Þ, with E the identity operation.Here, left-hand terms in the parentheses are point group operations and right-hand terms are translations."þ" and "−" refer to clockwise and counterclockwise rotations, respectively.Now considering only the point group operations, the MPG for the (001) surface is 3, which is indeed FM compatible.Specifically, a generically oriented magnetization pseudovector m ¼ ðm x ; m y ; m z Þ transforms in the following way under the point group operations in 3: Note that we have chosen our Cartesian coordinate system such that xk½110 and zk½001.From Eqs. ( 15)- (17), it is clear that the only invariant magnetization component is m z , which is parallel to the surface normal [001].Therefore, the magnetic dipole per unit area on the (001) surface must be oriented along [001].
[001] is known to be the groundstate polarization direction of the Néel vector [51][52][53], and from inspection of the electrostatically stable and polar surfaces in Figs.6(a) and 6(b) we see that the (001) surface always has an uncompensated ferromagnetic layer of Cr moments simply by terminating the ground-state magnetic order with vacuum.We can also check that this surface magnetization is roughness robust by going back to the R 30 c 0 bulk MSG and looking for operations that leave the direction of n unchanged but have a component of translation perpendicular to the surface.This is indeed the case for the three combined mirror plane-time-reversal operations for which n lies in the mirror plane: ðσ 0 ½100 j0; 0; 1 2 Þ, ðσ 0 ½110 j0; 0; 1 2 Þ, and ðσ 0 ½010 j0; 0; 1 2 Þ, where the prime denotes the additional action of Θ.Note that we use a convention where σ ½abc is the mirror plane perpendicular to the ½abc direction.Now we have to see how the magnetization behaves under these operations: Therefore, although operations ( 18)-( 20) take one (001) surface through a half-lattice-vector step to a symmetrically equivalent surface, all three operations leave m z , the FMcompatible direction of magnetization in the surface MPG, unchanged.Then, we also know that this magnetization is roughness robust.Note that this group-theory formalism makes no assumptions about the specific atomic termination and the corresponding value of (001) surface magnetization.The identification of this surface as uncompensated and roughness robust is termination independent, implying that surface magnetization cannot vanish for any (001) termination.Quantitative information about the value of surface magnetization for an experimentally realistic surface requires calculation of the multipolization array for chromia.
It turns out that, due to the alternating up-down-up-down order of the Cr moments, the only possible values of ME multipolization along this direction are those for the atomic terminations discussed in Sec.III A and shown in Figs.6(a) and 6(b) (11.9μB =nm 2 and −2.4μ B =nm 2 , respectively; for the opposite bulk AFM domain, the signs of both values are reversed).Thus, (001) Cr 2 O 3 has two symmetrically distinct terminations with unequal but finite magnitudes of surface magnetization, as required by the roughnessrobust identification from group theory.In practice, the termination corresponding to an electrostatically stable surface is energetically favorable, so experimentally all atomic steps are symmetrically equivalent to Fig. 6(a) with the corresponding surface magnetization.
Before closing this section, we mention one caveat regarding uncompensated, roughness-robust surface magnetization in ME AFMs.In spite of its vast technological potential and the fact that this category should exist for at least some surfaces of any linear ME AFM, to our knowledge the only roughness-robust, uncompensated surface magnetization which has been directly measured to date is our example case, (001) Cr 2 O 3 .The reason for the paucity of experimentally confirmed cases of uncompensated surface magnetization in ME AFMs remains unclear.A pragmatic possibility is that (001) Cr 2 O 3 has dominated interest and attention due to its relatively high Néel temperature compared to other MEs, making it and its doped variants promising candidates in applications [18,54,55].Nevertheless, given the myriad uses of such surface magnetization and the fact that the surface Cr ions on the (001) surface are quite weakly coupled to bulk near the Néel temperature [44], pursuing other ME AFMs as well as different surfaces in chromia could reveal new candidates with more robust uncompensated surface magnetization at room temperature.

Uncompensated, roughness-sensitive:
(111) NiO and (001) We turn next to surface magnetization that is magnetically uncompensated but which averages to zero in the presence of surface roughness (roughness sensitive).Let us again consider the minimal symmetry requirements of the bulk MSG in order for such surface magnetization to exist on any Miller plane.In contrast to uncompensated, roughnessrobust surface magnetization described in Sec.III B 1, uncompensated, roughness-sensitive surface magnetization corresponds to a centrosymmetric local-moment ME multipolization array.Without such a centrosymmetric multipolization array, there cannot be symmetry-equivalent surface steps with opposite directions and equal magnitudes of magnetization.
As we mention in Sec.III A, however, the multipolization array can take the centrosymmetric form of Eq. ( 12) for specific directions even if the magnetic lattice breaks I symmetry, provided there is another operation in the MSG that constrains the array to be centrosymmetric in the direction of interest.Thus, in principle, uncompensated roughness-sensitive surface magnetization can occur in AFMs with magnetic lattices having both broken and preserved I.We point out that, for all examples we investigate, the magnetic lattice has I symmetry.Indeed, the examples we work through explicitly (NiO and Fe 2 O 3 ) both have inversion symmetry in the bulk MSG.PbNiO 3 , a rhombohedral AFM insulator which has attracted interest due to possible ferroelectric properties [56,57], is an example of an AFM with broken inversion symmetry in the bulk MSG but a centrosymmetric lattice of Ni ions, which allows for roughness-sensitive, uncompensated surface magnetization along certain directions.The fact that we have not been able to identify any examples of roughness-sensitive surface magnetization in AFMs with I-broken magnetic lattices implies that the situation is likely somewhat contrived.
We now move on to the requirements for Θ in order to have uncompensated, roughness-sensitive surface magnetization.As discussed previously, Θ reverses the direction of any magnetization pseudovector m while leaving the polar vector corresponding to a surface normal n unchanged [that is, Eqs. ( 2) and ( 4) are fulfilled].Obviously then, if pure time reversal, or Θ combined with any translation parallel to the surface, is a symmetry of the bulk MSG, it will be preserved in the surface MSG obtained in Fig. 2(b).This means the surface MPG will not be FM compatible, and, hence, the Miller plane in question cannot host surface magnetization.
Recall, however, that, in identifying the surface MSG and corresponding MPG in Fig. 2, we exclude symmetries in the bulk MSG which include a component of translation perpendicular to the surface, parallel to n.Therefore, if a symmetry ðΘjtÞ in the bulk MSG has t • n ≠ 0 [the equivalent statement to Eq. ( 3)], it is possible for the surface of interest ðhklÞ⊥ n to have a FM-compatible MPG in the atomically smooth limit.Thus, AFMs which have ðΘjtÞ as a symmetry of the MSG can have uncompensated, roughnesssensitive surface magnetization on surfaces for which t • n ≠ 0 for all t which Θ is coupled to.This implies that type-IV AFMs, which by definition contain ðΘjtÞ as a symmetry [49], can still have roughness-sensitive surface magnetization.It is important to note that, because time reversal is preserved macroscopically in such AFMs, any net ME response is forbidden.Finally, we mention that type-IV AFMs are implicitly excluded in the formalism of Ref. [2].
We now give examples of particular surfaces in two bulk AFMs which have uncompensated, roughness-sensitive surface magnetization.We first focus on NiO.NiO is an insulating AFM with rocksalt structure.The nonmagnetic, face-centered cubic unit cell is shown in Fig. 9 We go through the same procedure as for (001) Cr 2 O 3 , where we first identify the space group symmetries in the bulk MSG which leave nk½111 unchanged except for possible translations parallel to the surface.C c 2=c has a monoclinic axis along ½1 10 and, thus, contains twofold rotations about ½1 10 and perpendicular mirror planes σ ð1 10Þ combined with several linearly independent translations.I and Θ combined with translations are also symmetries.The operations which leave n invariant modulo t⊥ n are ðEj0Þ, ðEj 1 2 ; 0; 1 2 Þ, ðσ 0 ½1 10 j0Þ, and ðσ 0 ½1 10 j 1 2 ; 0; 1 2 Þ [60].Now considering just the point group operation σ 0 ½1 10 , we can look at how magnetization transforms under this symmetry, taking xk½ 1 1 2, yk½1 10, and zk½111: Because xk½ 1 1 2, magnetization on the (111) surface is allowed along the direction of the bulk Néel vector, consistent with Fig. 9(b).Note that from Eq. ( 21) we see that an additional out-of-plane component of surface magnetization can develop along zk½111 via canting.Indeed, the surface MPG which corresponds to the identity plus a single σ 0 operation is m 0 , which is a FM-compatible MPG as expected.
We now go back to the bulk MSG and check for symmetries which reverse the surface magnetization while translating the polar vector n perpendicular to the (111) surface.This is the case for the symmetry ðΘj1; 1; 1Þ [61].The action of this symmetry is shown by the dashed black arrow in the structure on the left in Fig. 9(b), and it connects the left-hand surface with surface magnetization along ½11 2 to the symmetrically equivalent (111) surface on the right-hand side with surface magnetization along ½ 1 1 2. Evaluating MLM xz using Eq. ( 10) in a coordinate system with xk½ 1 1 2 and zk½111 and assuming a magnetic moment of 2μ B for Ni 2þ yields an uncompensated (111) surface magnetization of −13ðμ B =nm 2 Þ (þ13ðμ B =nm 2 Þ) along ½ 1 1 2 for the left (right) surface in Fig. 9(b).These two terminations are, in fact, related by a multipolization increment along the [111] direction, and it is not possible to create an atomically smooth (111) termination with a different magnitude of surface magnetization.This is the result of (111) NiO having a centrosymmetric halfincrement-containing multipolization array; see Eq. (II A).In the presence of surface roughness (again, assuming the idealized nonreconstructed surface), the two energetically degenerate surfaces in Fig. 9(b) occur with equal probability and the surface magnetization averaged to zero.
Roughness-sensitive, uncompensated surface magnetization is also possible for AFM planes corresponding to zero-containing multipolization arrays.We illustrate this for the case of (001) Fe 2 O 3 .Below the Néel temperature of approximately 955 K, the Fe moments lie in the x-y plane with a weak FM canting, but we focus on the lowtemperature phase of Fe 2 O 3 below 260 K where the moments are oriented along [001] [28].Fe 2 O 3 is isostructural with Cr 2 O 3 , but its "up-down-down-up" magnetic order along the [001] direction, shown in Fig. 9(c), does not break inversion symmetry, and so the linear ME effect is not allowed (although the symmetry of Fe 2 O 3 does allow for an antiferroically ordered local linear ME response [32]).Moreover, the full bulk MSG of lowtemperature Fe 2 O 3 , R 3c [167.103], has I symmetry.Therefore, any uncompensated surface magnetization must correspond to a centrosymmetric multipolization array.Searching for the operations in R 3c which leave nk½001 invariant modulo translations parallel to the surface, we get the same set of clockwise and counterclockwise threefold rotations [Eqs.( 16) and ( 17 However, if we now search for symmetries in R 3c which leave the direction of nk½001 invariant but translate n perpendicular to the (001) surface, we get mirror planes ðσ ½100 j0; 0; 1 2 Þ, ðσ ½110 j0; 0; 1 2 Þ, and ðσ ½010 j0; 0; 1 2 Þ.Note that these symmetries are analogous to those for (001) Cr 2 O 3 in Eqs. ( 18)- (20), except that in this case they do not contain time reversal.Therefore, for (001) Fe 2 O 3 these operations do switch the sign of the surface magnetization along zk½001, meaning that (001) Fe 2 O 3 surface magnetization is roughness sensitive, in contrast to roughness-robust surface magnetization in (001) ME Cr 2 O 3 .The mirror plane operations connecting surfaces with equal magnitude and opposite signs of magnetization are shown by the dashed black arrows for the two pairs of inequivalent terminations in Fig. 9(c).
If we consider the (001) Fe 2 O 3 surface magnetization from the perspective of its multipolization array, evaluation of MLM zz assuming a local moment of 5μ B for Fe 3þ yields three distinct values each differing by a multipolization increment: þ24μ B =nm 2 , 0μ B =nm 2 , and −24μ B =nm 2 .Further [001] translations of the topmost Fe moment on the (001) surface simply yield a repeat of these three values.Therefore, the multipolization array for (001) Fe 2 O 3 is zero-containing.From the left-hand pair of structures in Fig. 9(c), it is apparent that the finite AE24μ B =nm 2 surface magnetization corresponds to the polar, oxygen-terminated (001) surface.The surfaces with positive and negative magnetization are symmetrically equivalent due to the mirror plane combined with a [001] translation.On the right-hand side in Fig. 9(c), we see that the multipolization value of zero, in fact, corresponds to the electrostatically stable (001) termination.There are again two such surfaces which are symmetrically equivalent and connected by a mirror plane, but in this case the surface magnetization for both is zero.Therefore, in the case of a zero-containing centrosymmetric multipolization array, every inequivalent termination contains a pair of degenerate surfaces with equal magnitude but opposite sign of surface magnetization.Since degenerate surfaces occur with equal probability, any finite contributions would cancel.Experimentally, of course, only the electrostatically stable surface on the right-hand side is likely to occur, corresponding to the zero value of the multipolization array in this case.
In part II, we have built on our work in Ref. [28], which proposed the local-moment ME multipolization as a bulk descriptor of uncompensated surface magnetization.We have further shown that the symmetry of the multipolization array (centrosymmetric or noncentrosymmetric) immediately identifies whether the surface magnetization is roughness robust or roughness sensitive and additionally yields a prediction of the quantitative values of surface magnetization, in the absence of reconstruction, in contrast to the purely qualitative group-theory procedure in part I. Next, we move on in part III to characterize the induced forms of surface magnetization using symmetry arguments combined with ab initio calculations.

IV. PART III: INDUCED SURFACE MAGNETIZATION IN TERMS OF ATOMIC-SITE MULTIPOLIZATION
In this part, we explore the connection between surface magnetization and the bulk ME response in detail for nominally compensated surfaces.In these cases, the induced surface magnetization is the dominant effect, since there are no uncompensated magnetic dipoles in the surface plane.We start by reviewing the bulk ME responses at first, second, and third order in the electric field, and outlining how these can be conveniently described using the multipoles of the magnetization density introduced in Sec.III A. Next, we discuss as concrete examples ð 120Þ and (100) Cr 2 O 3 and (110) FeF 2 , showing how to link the allowed bulk ME and magnetic multipoles to the induced surface magnetism.

A. Prelude: Local ME responses and atomic-site multipoles
As introduced in the introduction and Sec.III A, by definition ME materials show a net change in magnetization, δM, when an external electric field E is applied or, vice versa, exhibit a net change in polarization, δP, in the presence of an external magnetic field H.The lowest-order, most well-known ME response is the linear ME effect [17], whereby the net change in magnetization is linear in the electric field's strength; see Eq. (8).
The linear ME effect can be conveniently recast in terms of microscopic indicators, called ME multipoles [25], corresponding to the entries of the ME multipole tensor M ij [26] introduced in Sec.III A, Eq. ( 6).ME multipoles have a one-to-one link to the linear ME tensor, with a nonvanishing ij entry of the ME multipole tensor M implying an ij ME response.
As discussed in Sec.III A, M ij can be decomposed into a sum of two terms.One is an origin-dependent, multivalued contribution, referred to as the "local-moment term" M LM ij , and captures the inversion-breaking asymmetries arising from the arrangement of the atomic magnetic dipoles.We used this first term in part II as a bulk indicator of surface magnetism for magnetically uncompensated surfaces.The other part of M ij is the atomic-site contributions, which describe the inversion-breaking asymmetries in the local magnetization density around the individual ions.In the following, we consider this atomic-site component, which we recently showed to be useful in predicting the local linear ME effect [32].Here, by local linear ME effect, we denote the change δm in the local atomic magnetic moment, induced by an external electric field: A nonzero local ME multipole implies a nonvanishing corresponding local linear ME response α loc ij .Whether this results in a net linear ME response α ij depends on the arrangement of the local ME multipoles in the unit cell: Specifically, ferroically ordered ME multipoles entail a net linear ME response, whereas antiferroically ordered ME multipoles imply a linear anti-ME response and, in turn, a vanishing net α ij .As a consequence of the previously mentioned symmetry correspondence between surface magnetization and ME responses, we soon see that a ferroic arrangement of these atomic-site multipoles in the surface plane indicates the presence of induced surface magnetization.On the other hand, whether the sign of in-plane ferroic order of the atomic-site multipoles switches (does not switch) between adjacent planes along the direction perpendicular to the surface distinguishes roughness-sensitive (roughness-robust) induced surface magnetization.
We also find that atomic-site multipoles of higher order than those associated with the linear ME effect can be important in explaining surface magnetism.In particular, we make the connection to the second-order and third-order ME responses, where the induced magnetization is bilinear and trilinear, respectively, in the applied electric field: Similarly to the linear ME response, the second-order response, described by the β ijk tensor, has been recently recast in terms of higher-order ME multipoles of the magnetization density, specifically magnetic octupoles [22]; see Eq. (7).In Appendix E, we report the irreducible spherical components of the magnetic octupole tensor O ijk , which is used later in the discussion.Likewise, the third-order ME response is associated with inversionbreaking magnetic hexadecapoles, although a detailed analysis of the connection, to our knowledge, is still missing in the literature.

B. Induced, roughness-robust surface magnetization
We move now to cases of magnetically compensated planes in the bulk AFM ground state for which roughnessrobust magnetization is induced at the surface.Again, let us consider the minimal symmetry requirements for this scenario in terms of I and Θ symmetries.As discussed in Secs.III B 1 and III B 2, roughness-robust surface magnetization cannot exist for any Miller plane when ðΘj0Þ or ðΘjtÞ, where t is any translation, is a symmetry of the bulk MSG.Therefore, we know right away that we must have broken Θ in the bulk.
We next consider inversion symmetry I.In Secs.III A and III B 1, we show that broken I in the bulk MSG is necessary to obtain the noncentrosymmetric form of the local-moment multipolization array in Eq. ( 14) necessary for uncompensated, roughness-robust surface magnetization.However, for surfaces which are magnetically compensated, the local-moment contribution [Eq.(10)] for the corresponding multipolization component is always zero.Therefore, for compensated surfaces, there is no contribution to the surface magnetization from the local-moment ME multipolization array.The lack of a local-moment multipolization array contribution for compensated surfaces with induced magnetization implies that I may be either present or absent in the bulk MSG [62].We show this explicitly by leveraging the group-theory analysis, combined with first-principles DFT calculations, for three example AFM surfaces, two with broken inversion in the bulk and one with a centrosymmetric bulk MSG.In addition, we show that, in the roughness-robust case, induced surface magnetization is associated with atomicsite multipoles of the magnetization.Specifically, these are the inversion-broken atomic-site ME multipoles and hexadecapoles for ð 120Þ and (100) Cr 2 O 3 and inversionsymmetric atomic-site magnetic octupoles for the (110) surface of centrosymmetric FeF The surface normal n is parallel to [010], which is a twofold rotational axis in the bulk MSG.If we go through all space group operations in R 30 c 0 [see Fig. 10(b)] and select those which leave nk½010 invariant modulo translations ⊥ n, we find (besides the identity) only the twofold rotation about the surface normal, ð2 ½010 j0; 0; 1 2 Þ.Adopting a rotated basis where xk½010 and zk½001, the magnetization transforms in the following way: implying that a finite surface magnetization is allowed along x.Moreover, there are no symmetries in R 30 c 0 which flip the sign of m x and translate n perpendicular to the surface, so the surface magnetization is roughness robust.Note that the only direction which is FM compatible for the ð 120Þ surface, xk½010, is not parallel to the bulk [001] Néel vector polarization but, rather, along the surface normal.Therefore, surface magnetization is allowed by symmetry to develop only via canting along [010].Interpretation in terms of linear ME effect.Next, we rederive from a different standpoint the results obtained in part I using the group-theory-based prescription to identify surface magnetism.In particular, motivated by the equivalence between surface magnetism and the bulk ME response, we address the magnetization induced in the bulk by an electric field parallel to the surface normal.Before doing so, we remind the reader that bulk Cr 2 O 3 allows for a net linear ME response identified by the ME tensor which is proportional to the ME multipole tensor.Since we are interested in the local atomic response to an electric field, here we consider the local linear ME response for the Cr atoms, which reads Compared to the bulk response [Eq.( 26)], the local ME tensor shows additional off-diagonal xy entries, since the site symmetry of the Cr Wyckoff position ( 30) is lower than the point group symmetry ( 30 m 0 ).We use overlined symbols in Eq. ( 27) to indicate antiferroically ordered local ME multipoles (and the components of the ME response associated with them).Now we consider the local response to an electric field Ek½010.The [010] direction makes a 120°a ngle with the x axis in the hexagonal setting [Fig.10(a)]; thus, we conveniently rotate the reference framework to align the [010] axis along the x axis.In this way, the ME response along [010] to an electric field parallel to [010] corresponds to the xx entry of the rotated tensor.The rotation transformation for α loc reads where the prime superscript refers to the rotated framework and is the matrix representing a clockwise rotation of angle θ ¼ 120°.Equation ( 28) implies α 0loc xx ¼ α loc xx ; thus, the magnetic moment induced along x 0 k½010 is δm x 0 ¼ α loc xx E. This result is formally equivalent to the surface magnetization predicted by group-theory arguments earlier in this section.
(100) surface, symmetry analysis.The (100) surface of chromia, similarly to the ( 120) surface, is obtained by cutting bulk Cr 2 O 3 along a plane perpendicular to the basal hexagonal plane.Its normal nk½210 is related to the normal of ( 120), that is, [010], by a clockwise 90°rotation about the z axis; see Fig. 10(a).Following the group-theory-based guidelines to identify whether this surface allows for roughness-robust surface magnetism, we first identify which symmetry operations among those of the magnetic space group R 30 c 0 of chromia [Fig.10(b)] leave n invariant up to a fractional translation perpendicular to n itself.By inspection, one finds that the operations fulfilling this requirement are the identity E and the mirror plane perpendicular to [010] combined with time reversal, σ 0 ½010 .After defining a rotated reference framework (x, y, z), with xk n, yk½010, and zkz, the action of σ 0 ½010 on a generic magnetic moment m ¼ ðm x ; m y ; m z Þ is which leaves m x and m z invariant.Thus, (100) Cr 2 O 3 is compatible with ferromagnetism both along the surface normal and along the [001] direction, within the surface plane.This is distinct from the ( 120) surface analyzed earlier, which was compatible with ferromagnetism only along the surface normal.
Interpretation in terms of linear ME response.Now, we interpret the surface magnetism in terms of the bulk linear ME response to an electric field along the (100) surface normal [210].Following the argument discussed above for the ( 120) surface, we take the local linear ME tensor α loc and we apply a clockwise rotation of angle θ ¼ 30°in order to align the [210] direction to x.We find α 0loc ¼ α loc , which implies After summing over the Cr atoms on the topmost surface atomic layer, δm y vanishes, since ᾱloc xy is antiferroically ordered, leaving only a nonvanishing net δm x .Thus, using the bulk linear ME response we are able to explain only the component of the surface magnetization perpendicular to the surface, and we are not able to predict the component parallel to the surface, along [001].This feature comes from the local ME response being isotropic in plane, which yields a contribution of surface magnetization perpendicular to the surface for any ðhk0Þ surface, perpendicular to the basal plane.In order to explain components of the magnetization parallel to the surface, we need to go beyond the linear ME response and account for higher-order, anisotropic responses, which we discuss next.
Higher-order ME-induced contributions.As discussed in previous sections, based on symmetry arguments we predict that the ( 120) and (100) surfaces of Cr 2 O 3 are not equivalent concerning the allowed components of net surface magnetization: Specifically, the ( 120) surface is compatible with ferromagnetism only along [010], perpendicular to the surface, whereas the (100) surface allows for ferromagnetism both along the [210] direction, perpendicular to the surface, and along the [001] direction, parallel to the Néel vector.As examined earlier in Secs.IV B 1 a and IV B 1 c, the isotropic linear ME response explains the component of surface magnetization perpendicular to the surface but does not explain the component parallel to zk½001, since a zx linear ME response is forbidden by symmetry in the bulk.In order to understand, in terms of ME responses, why such a component of the surface magnetization parallel to [001] exists and why it is allowed only in the family of (100) and equivalent surfaces, we need to analyze higher-order ME responses, beyond the linear one.For clarity, in the following, we adopt the right-handed coordinate system with xk½100, yk½120, and zk½001, shown in Fig. 10(a).As such, x lies on the twofold rotation axis 2 ½100 of the 3m point group, whereas y is on the mirror plane σ ½100 , and z is the 3 roto-inversion axis [see Fig. 10(b)].
We start by analyzing the second-order ME response which, as mentioned earlier in Sec.IVA, is conveniently recast in terms of magnetic octupoles.In Cr 2 O 3 , the octupoles allowed by symmetry are O 3 , Q ðτÞ z 2 , which follow a (þ þ −−) order [here, the signs refer to the octupole component on Cr 1 , Cr 2 , Cr 3 , and Cr 4 atoms as labeled in Fig. 8(a)], and O −3 , O 0 , and t ðτÞ z , which follow a (þ − þ−) order.Note that here we follow the same naming convention for magnetic octupoles as in Ref. [22].Both families of octupoles order antiferroically; hence, despite contributing to the local second-order ME response, they do not provide any net response.The implication for surface magnetism is that an induced local atomic magnetic moment, bilinear in the effective electric field generated by the surface termination, is expected for the Cr atoms at the surface.However, for (100) and ð 120Þ surfaces these induced atomic magnetic moments arrange antiferroically, since the magnetic octupoles of the four Cr atoms appearing at the surface plane are antiferroically ordered and, thus, do not yield any net surface magnetization.
Next, we consider the third-order ME response.To simplify the discussion, we do not describe in detail the induced local atomic response but limit ourselves to the net third-order ME response, for which the induced net magnetization δM is trilinear in the external electric field E: with γ ijkl the rank-4 third-order ME response tensor.Note that, as mentioned earlier, γ ijkl is associated with magnetic hexadecapoles, defined as H ijkl ¼ R Ω r i r j r k μ l ðrÞd 3 r, but, since a detailed analysis of such connection is so far missing in the literature, here we refer simply to the third-order ME response instead.
Since γ ijkl breaks inversion symmetry, such a response is symmetry allowed in noncentrosymmetric Cr 2 O 3 .As we want to explain the surface magnetization along the z (k½001) direction in the (100) and symmetry-equivalent surfaces, we focus on the net magnetization induced along z by a general electric field in the basal plane; thus, we consider γ zjkl , with j ¼ fx; yg, k ¼ fx; yg, and l ¼ fx; yg.We use the MTENSOR utility of the Bilbao Crystallographic Server [63][64][65], to write γ zjkl in terms of its independent parameters as First, we analyze the ð 120Þ surface.We note that x (k½100) is symmetry equivalent to the [010] direction normal to ð 120Þ [see Fig. 10(a)], since the two directions are connected by a 120°rotation; thus, the net magnetization along z induced by an electric field along x is equivalent to the surface magnetization along z for the ð 120Þ surface.Since γ zxxx ¼ 0, our tensor analysis predicts no surface M z for this surface, in agreement with what we obtained earlier in Sec.IV B 1 a from group-theory arguments.Next, we consider any (hk0) surface, with normal vector n lying in the hexagonal basal plane [see Fig. 10(a)].For convenience, we apply a rotation to the tensor γ, such that it aligns the electric field Ek n, perpendicular to the surface, with the x axis, similarly to our earlier procedure for the linear ME contribution.To keep the discussion general, we apply a clockwise rotation R about the [001] axis, for a generic angle θ: The tensor in the rotated framework is computed as follows: As we are interested in the zxxx response in the rotated framework, we set i 0 ¼ z and j 0 ¼ k 0 ¼ l 0 ¼ x.From Eqs. ( 33)-( 35), we obtain We note that, in contrast to the linear ME response, the third-order response is anisotropic in plane.We show γ 0 zxxx in Fig. 12, where we arbitrarily set c zxxy to 1.We note that γ 0 zxxx as a function of θ is periodic with period 120°, arising from the bulk threefold rotational symmetry.In addition, importantly, γ 0 zxxx ≠ 0 for the (100), (010), and ( 110) surfaces, which correspond to θ ¼ 30°; 90°; 150°, respectively, whereas γ 0 zxxx ¼ 0 for the (2 10), (110), and ( 120) surfaces, corresponding to θ ¼ 0°; 60°; 120°, respectively.Thus, based on the analogy between bulk ME responses and surface magnetism, the anisotropic third-order ME response explains the behavior of the surface magnetization along z that we predicted earlier using group theory for the (100) and ( 120) surfaces.
Ab initio calculations.So far, we have used symmetry arguments to show that surface magnetization can be induced even when the magnetic dipoles in the plane are magnetically compensated.Now, we corroborate our predictions using ab initio calculations.This subsection has two purposes: first, to demonstrate with quantitative DFT calculations that a finite magnetization, obtained via canting of the magnetic moments, on the ð 120Þ and (100) surfaces of Cr 2 O 3 lowers the total energy compared to the surface with zero magnetization and, second, to strengthen the link between surface magnetization and bulk ME effects by testing with first-principles calculations whether the sign change of the ME response upon bulk ME domain reversal, a core ingredient of ME domain selection using ME annealing [2,55], has its counterpart in surface magnetism.If the magnetization induced on a surface can be indeed interpreted as a consequence of a ME effect from an effective electric field due to the surface termination, we expect the surface magnetization to reverse in the opposite bulk magnetic domain.To test this assumption, we consider both ( 120) and (100) surfaces obtained from both bulk magnetic domains, as shown in Fig. 13.We refer to these domains as the in-pointing and the out-pointing domains, depending on whether the magnetic moments of each pair of nearest-neighbor Cr atoms point toward each other ("in") or away from each other ("out").
In our DFT calculations, whose technical details are described in Appendix A, we fix the Cr moments in the center two (four) layers of our four-layer (six-layer) ð 120Þ-oriented [( 100  The change in total energy per formula unit with respect to the energy with 0°canting (no out-of-plane surface magnetization) and the net out-of-plane surface magnetization, as a function of canting angle, for both domains, are shown in the bottom in Figs.14(a) and 14(b), for the ð 120Þ and (100) surfaces, respectively.Beginning with the ð 120Þ case [Fig.14(a)], first we notice that, rather than an energy minimum for zero canting, the slab energy instead reaches a minimum when the surface moments are canted around 0.25°with respect to [001].While the energy difference is tiny (< 1 μeV), the penalty energy for constraining the moments (defined in Appendix A), plotted with red stars and green circles for the in-pointing and out-pointing domains, respectively, is orders of magnitude smaller.This, combined with our stringent convergence criterion (we stop the electronic minimization when energy differences between successive iterations differ by 10 −8 eV or less), makes us confident in the physical significance of the energy lowering.For 0.25°c anting, we obtain a ð 120Þ surface magnetization with magnitude ∼0.12μ B =nm 2 , calculated by taking the projections of all 12 Cr surface moments onto the [010] normal direction and dividing by the ð 120Þ surface area.Second, we remark that the energy curves corresponding to the in-pointing and the out-pointing domains (red and green stars) are the mirror images of one another with respect to zero canting angle; hence, their minima correspond to oppositely oriented surface magnetizations.This implies that the surface magnetization reverses if we switch from one domain to the other, as we expect from the interpretation in terms of ME responses.
The results for the (100) surface, shown in Fig. 14(b), are qualitatively similar to those for ð 120Þ but quantitatively slightly different, as the canting perpendicular to the surface is predicted to be approximately 0.5°.Although the induced moment per Cr correspondingly doubles for the (100) surface compared to that for ( 120), the density of surface Cr per unit area decreases, such that the out-ofplane surface magnetizations in units of μ B =nm 2 for the ð 120Þ and (100) surfaces are almost identical, as can be seen in Fig. 14.Again, the surface magnetizations for the two opposite domains on the (100) surfaces are oppositely aligned with respect to each other, as expected based on the connection to the ME response.
We now comment on the additional surface magnetization along [001] which we expect for the (100) surface based on our symmetry analysis in Secs.IV B 1 c and IV B 1 e.In our DFT calculations, we impose a spin canting perpendicular to the surface, and we constrain only the direction of the atomic magnetic dipoles to be in the plane spanned by the surface normal and the [001] directions, but we do not constrain their magnitude.In this way, we, in principle, allow the magnetic moments to arrange in a way that develops a net component along [001] when summed over the Cr atoms.As discussed in Sec.IV B 1 a, for the case of the ð 120Þ surface, finite magnetization along [001] is forbidden by the symmetry of the MPG.Thus, the sublattice projections along [001] remain equivalent, as indicated by the same-color moments for the slab cartoon in Fig. 14(a).The only surface magnetization induced is that along the ½ 120 surface normal, as plotted in the bottom in Fig. 14(a).On the other hand, based on the (100) surface MPG, we expect an additional component of surface magnetization parallel to the Néel vector along [001].Indeed, this is confirmed by our DFT calculations, in which the Cr magnetic moments oppositely aligned along [001] on the top and bottom layers of the surface become inequivalent and show a difference in magnitude of approximately 0.03μ B (depicted pictorially by the differentcolored surface sublattices) in Fig. 14(b).The sign of the surface magnetization along [001] also switches with opposite domains, as expected from analogy to the thirdorder ME response.The subsequent in-plane induced surface magnetization is about half the size of the outof-plane surface magnetization corresponding to the energetic minima in Fig. 14(b).
Summary.Before ending this section, we comment briefly on the practical relevance of our findings for ð 120Þ and (100) Cr 2 O 3 .Admittedly, the value of induced surface magnetization, 0.12μ B =nm 2 , is small.We stress, however, that our calculations are not a rigorous indicator of the quantitative, experimental induced magnetization that one would expect.Our goal was a proof-of-principle demonstration of energy lowering for a finite surface magnetization rather than determination of the groundstate magnetization for a ð 120Þ or (100) surface.For thicker slabs with a larger number of central bulklike layers, it is likely that a surface could be stabilized further by canting the moments for several of the outermost layers rather than just a single layer, thus modifying the total magnetization compared to the values we obtain here.Additionally, for exchange bias applications where the total surface magnetization summed over all unit cells at the FM-AFM interface is more relevant than the surface magnetization per unit cell, the small induced magnetization on ð 120Þ and (100) Cr 2 O 3 could still have substantial effects.Encouragingly, recent independent experiments using magnetic circular dichroism and magnetotransport have actually detected a nonzero surface magnetization on the nominally compensated surfaces of Cr 2 O 3 perpendicular to (001) [14,66], suggesting that our theoretical prediction of induced surface magnetization in this case can be be experimentally detected.
Moreover, while the magnitude is small compared to the uncompensated case in Sec.III B 1, we still expect this surface magnetization to be much larger than the bulk magnetization induced by the transverse ME effect (E⊥½001) for realistic electric field strengths.Previous first-principles calculations of the lattice-mediated contribution to the transverse ME response indicate that, for a typical electric field strength of 1 V=nm, the induced magnetization in Cr 2 O 3 is about 0.002μ B per Cr along E [32,67], corresponding to an approximately 0.04°canting angle.In contrast, the surface magnetization toward the surface normal which we find for the ð 120Þ and (100) surfaces of chromia for canting angles of 0.25°and 0.5°t ranslates to an induced 0.011μ B and 0.022μ B , respectively, per Cr ion.Thus, the surface magnetization per Cr induced by a vacuum-terminated surface is substantially larger than that induced due to the bulk ME effect in Cr 2 O 3 for a reasonably sized electric field.This is not surprising when considering that the effective electric field of a vacuum-terminated boundary should correspond roughly to the scale of the material crystal field, which can be orders of magnitude larger than an experimentally achievable electric field [2].

Induced surface magnetization in (110) surface of FeF 2
We next consider induced, roughness-robust surface magnetization in a centrosymmetric AFM.We choose FeF 2 , a rutile-structure insulating AFM with a Néel temperature of approximately 80 K [68], for a few reasons.First, induced magnetization on the (110) FeF 2 surface was recently detected experimentally by analyzing the giant magnetoresistance at the interface of an FeF 2 -Cu-Co spin valve [13].Second, as mentioned in the introduction, the (110) surface of FeF 2 as well as isostructural (110) MnF 2 exhibit large exchange bias experimentally, despite the fact that these surfaces are magnetically compensated when considering the bulk AFM ground order [8,11,19,69].Many possible explanations have been proposed, including structural matching between the AFM and FM surfaces, noncollinear coupling at the interface, alignment of interface moments during field cooling, or piezomagnetismgenerated uncompensated interface moments [8,11,19].However, given the recent theoretical and experimental evidence that these surfaces have a symmetry-based equilibrium magnetization, a connection between induced surface magnetization and the exchange bias for (110) FeF 2 seems likely.
FeF 2 crystallizes in the centrosymmetric tetragonal rutile structure described by the P4 2 =mnm space group (4=mmm point group).Its unit cell, shown in Fig. 15(a), contains two formula units.The two Fe atoms sit at (0, 0, 0) and (1=2, 1=2, 1=2) and are octahedrally coordinated by the surrounding F atoms.The FeF 6 octahedra are elongated along one axis [see Fig. 15(b)], and, importantly, for the two Fe atoms in the unit cell they are rotated by 90°about the [001] direction with respect to each other.Hence, despite belonging to the same Wyckoff orbit (2a) and being thus symmetry equivalent, the two Fe atoms are not connected by any pure translation.This has an important consequence when the magnetic ordering is taken into account, since the inequivalent fluorine environment of the Fe atoms makes the AFM order of FeF 2 , described by the magnetic space group P4 0 2 =mnm 0 , break time-reversal symmetry.
(110) surface, symmetry analysis.The bulk unit cell defining the electrostatically stable (110) surface consists of two Fe layers; see Fig. 16(a).Figure 16(b) shows the eight-layer slab used in our DFT calculations.Each layer stacked along the [110] direction contains one Fe ion per unit cell from the corner of the bulk unit cell and one from the center sublattice, such that the (110) surface is magnetically compensated with oppositely pointed in-plane moments.We now find the symmetries of the bulk MSG which leave nk½110 invariant modulo translations parallel to the surface.Such symmetry analysis for the existence of (110) surface magnetization was already discussed in Ref. [13], but we repeat it here briefly.In addition to the trivial identity operation, these are time reversal combined with a twofold rotation ð2 0 ½110 j0Þ about the ½110 surface normal, as well as two mirror planes ðσ ½001 j0Þ and ðσ 0 ½1 10 j0Þ.This corresponds to the FM-compatible surface MPG 2 0 mm 0 .We now check how magnetization transforms under these operations, taking a rotated Cartesian basis with xk½110, yk½1 10, and zk½001: Equations ( 37)-( 39) indicate that the only direction along which a magnetization can develop is zk½001.Thus, in contrast to ð 120Þ and (100) Cr 2 O 3 , for which out-of-plane surface magnetization developed via a relativistic canting of the surface moments perpendicular to the bulk Néel vector direction, for (110) FeF 2 the surface magnetization develops only parallel to the bulk Néel vector.
In order for a finite magnetization to develop parallel to the Néel vector, the moment magnitudes of the two sublattices in the surface layer must necessarily become inequivalent.This is allowed by symmetry, since none of the operations (37)- (39) in the surface MPG connect the corner sublattice moment to the center sublattice moment for any (110) layer.Moreover, as pointed out in Ref. [13], we can understand intuitively why the symmetry-allowed imbalance in sublattice magnitude might be particularly enhanced for the vacuum-terminated surface layers [layers 1 and 8 in Fig. 16(b)] by considering the F coordination of the surface Fe ions.This is more easily seen by rotating the (110) slab such that the [001] Néel vector polarization points into the page, as done in Fig. 16(a).From this, it is apparent that, for the electrostatically stable (110) surface, the upper F atom for the "corner" Fe sublattice is cut off (indicated by its greater transparency in the figure), whereas the "center" Fe sublattice retains its full bulk coordination of F ions.Thus, the chemical environments around the Fe ions on the surface layers are strongly dissimilar, implying that the inequivalence in the magnitude of their moments will be large for the surface layers.We note that, in contrast to the ð 120Þ and (100) surfaces of Cr 2 O 3 , a surface magnetization in (110) FeF 2 does not arise because of spin canting; thus, the mechanism is nonrelativistic in nature.From a symmetry point of view, the surface magnetization in (110) FeF 2 arises because the equivalence of the Fe sublattices is lifted by the symmetry lowering due to the surface effective electric field.We remark that the cases of surface magnetism arising from spin canting and from sublattice-equivalence lifting can be distinguished also with group-theory-based arguments, as explained in Ref. [2].
Interpretation in terms of ME effects.Next, we show how the surface magnetism for the (110) surface of FeF 2 is interpreted in terms of ME effects.Since FeF 2 is centrosymmetric, with the Fe atoms at inversion centers, a linear ME response, both local and net, is forbidden by symmetry.Instead, the lowest-order ME response is the second-order ME effect, in which the induced magnetization is bilinear in the applied electric field, M ∝ βE ⊗ E; see Eq. ( 23).In FeF 2 , this manifests both as a local, atomic response and as a bulk net effect due to the breaking of the time-reversal symmetry.Following Ref. [22], the local second-order ME response can be conveniently interpreted in terms of magnetic octupoles.(note that here we follow the naming convention of Ref. [22]).Overall, the full magnetic octupole tensor (see Appendix E) at individual Fe sites in FeF 2 reads We note that, besides being symmetric upon exchange of j and k by construction, O xjk and O yjk are equivalent on exchange of x and y, due to the tetragonal symmetry (with fourfold rotation axis parallel to z) of both the crystal and the magnetic structure.For the same reason, O zjk is not equivalent to O xjk and O yjk .
As a next step, we consider the local second-order ME response.The Fe magnetic moments δm induced by an external electric field E are given by The tensor β loc ijk describes the local second-order ME response and has the same symmetry properties as the magnetic octupole tensor [22]; thus, it can be written in terms of independent parameters in the following way: where the overlined (nonoverlined) symbols identify the antiferroic (ferroic) responses, i.e., those that have opposite (same) sign on the two Fe atoms.Following again the analogy between surface magnetism and ME effects, in order to study the magnetism of the (110) surface, we analyze the bulk ME response to an electric field along the [110] direction.Previously, we showed that this can be done by applying a rotation to the tensor, such that it brings the electric field from the [110] direction to the x direction.Thus, we apply a clockwise rotation R of angle θ ¼ 45°a bout the z axis: to the rank-3 tensor β loc , which accordingly transforms as follows: In the rotated framework, we are interested in the atomic magnetic moments induced by an electric field along x.From Eqs. ( 41)-( 44), we obtain which implies that the only second-order ME response allowed in this case is off diagonal and that it does not induce any spin canting but rather a change in size of the zero-field magnetic moment along the Néel direction zk½001.Interestingly, δm z contains both a contribution from the antiferroically ordered O 0 and t ðτÞ z and the ferroically ordered O −2 and Q ðτÞ x 2 −y 2 .As a consequence, the local magnetic moments induced on the two Fe atoms, are different in size.Thus, the second-order zxx ME response in this case is a ferri-ME response.We remark that the different response of the two Fe atoms stems from their geometrically inequivalent fluorine environment, which is also crucial for the time-reversal symmetry breaking and, consequently, of the nonrelativistic spin splitting of the electronic energy bands, similarly to MnF 2 [39].Since the (110) surface magnetization can be interpreted in terms of a second-order bulk ME response to an electric field perpendicular to the surface, our results suggest a net magnetic moment parallel to the (110) surface, which is due to inequivalent antiparallel magnetic moments, with different magnitude on the two Fe atoms, in agreement with the group-theory-based analysis discussed earlier.
Ab initio calculations.Next, we confirm our group-theory and second-order ME tensor analyses with first-principles DFT calculations, the details of which are described in Appendix C. First, we confirm that a finite surface magnetization does indeed develop along the [001] inplane direction for the (110) surface.We plot our calculated, site-projected Fe magnetization (in μ B ) as a function of layer number for corner and center Fe sublattices on the top and center subplots, respectively, in Fig. 16(c).The corresponding "layer magnetization" in μ B =nm 2 , calculated by summing the two sublattice moments in each layer and dividing by the cross-sectional (110) surface area, is plotted in the bottom in Fig. 16(c).For central layers 3-6, the magnitudes of the oppositely pointed moments are nearly constant and identical in size (3.12-3.13μB ), and, correspondingly, the bulk layers are almost perfectly compensated as seen in the bottom plot.For surface layers 1 and 8, however, we see an appreciable 0.008μ B difference in sublattice magnitudes, corresponding to a negative surface magnetization of −0.035μ B =nm 2 along [001].For our choice of AFM domain [Fig.15(a)], this corresponds to a larger moment magnitude on the center Fe sublattice.
Clearly, the precise numerical values for all three subplots in Fig. 16(c) are dependent on the specifics of the DFT parameters, such as the Hubbard U and Hund J as well as the radius of the integration sphere used to calculate the site-projected magnetization (all of which we report in Appendix C).Nevertheless, our results provide an unambiguous, quantitative confirmation of the symmetry analysis.Second, we confirm our predictions for the local and net bulk ME response with additional ab initio calculations on bulk FeF 2 .These are carried out by freezing in infraredactive phonon modes polarized along the [110] direction, corresponding to the atomic displacements induced by an electric field pointing along [110].For more details on computational methods and parameters, see Appendix C. According to our calculations, FeF 2 shows seven infraredactive phonon modes.Six of these transform as the E u irreducible representation of the 4=mmm point group: They are polarized in the ðx; yÞ plane, identified by the [100] and [010] directions (see Fig. 15) and are grouped into three families of degenerate modes.The remaining infraredactive mode transforms as the A 2u irreducible representation of 4=mmm; thus, it is polarized along the z direction.Since we are interested in the response to an electric field along the [110] direction, here we consider the contribution of the E u modes to the ME response.As a proof of concept, we focus on the contribution of the highest-frequency E u mode.We construct an appropriate linear combination of the two degenerate modes such that the resulting mode is polarized along [110].Subsequently, we freeze this mode into the structure and study the induced magnetic moments on the Fe atoms as a function of the mode's amplitude.Figure 17 summarizes our results.Black and red data points show the induced magnetic moment along z for the two Fe atoms.The response is quadratic in the mode amplitude and has a different magnitude for the two Fe atoms, thus rendering them inequivalent, as predicted by Eq. (46).FIG. 17. Change in the magnetic moment along z on the two Fe atoms in FeF 2 (red squares and black circles) and sum of the Fe magnetic moments along z (black triangles), induced by an infrared-active phonon mode polarized along [110].
The difference in responses on the Fe sites results in a net induced magnetic moment when summing over the two ions, as shown by the blue data points in Fig. 17, in agreement with our tensor and symmetry analyses above.
Summary.Before concluding our discussion of roughnessrobust, induced surface magnetization, we point out one additional exciting implication of surface magnetization for AFMs such as FeF 2 , in addition to its likely connection to exchange bias.A recently identified class of AFMs, namely, altermagnets, exhibiting nonrelativistic spin-split band structures along particular directions in reciprocal space [39,49,50,70] has generated great excitement due to their potential use in applications such as spin-charge conversion.Altermagnets require broken ΘI symmetry in the bulk MSG, and, thus, linear ME AFMs with preserved ΘI, such as Cr 2 O 3 , are not altermagnets [71], whereas centrosymmetric FeF 2 and its isostructural counterpart MnF 2 [49,50] are.
Crucially for us, the directions in the Brillouin zone along which nonrelativistic spin splitting is allowed correspond to k points with a little group of k (the subgroup of the MPG leaving k invariant modulo a reciprocal vector) that does not reverse the spin direction [50].Note that this is similar to the definition of a FM-compatible surface MPG, with the surface normal n corresponding to the wave vector k.Of course, there is not a one-to-one correspondence, first because, unlike n, k switches sign with time reversal, and, second, because real and reciprocal space directions are not parallel for nonorthogonal lattices.For cubic and tetragonal AFMs such as FeF 2 , however, the paths in reciprocal space with nonrelativistic spin splitting coincide exactly with surfaces allowing induced surface magnetization.Specifically, in FeF 2 and MnF 2 , the high-symmetry paths with nonrelativistic spin splitting are Γ-M and Z-A, corresponding to k ¼ ðu; u; 0Þ and k ¼ ðu; u; 1=2Þ, as well as parallel paths at other constant k z values.These directions are parallel to the real space (110) surface normal for which we and previous authors have demonstrated the existence of roughness-robust induced surface magnetization.Therefore, experimental detection of surface magnetization could be an alternative indicator of nonrelativistic spin splitting, and vice versa.

C. Induced, roughness-sensitive surface magnetization
We note that it is also possible to have AFM surfaces that have an induced magnetization which averages to zero in the presence of roughness.Given the combination of small surface magnetization in the induced case and roughness sensitivity, this fourth category of surface magnetization is likely not very useful.Moreover, thus far, we have not found an explicit material example of this category.We nevertheless discuss the symmetry and multipole requirements for roughness-sensitive, induced surface magnetization and leave the analysis of a concrete example for future work.
First, let us note the minimal bulk symmetry requirements in terms of I and Θ which allow for induced, roughness-sensitive surface magnetization.Based on analogous arguments to those given in Sec.IV B, this surface magnetization can occur in AFMs with either broken or preserved inversion symmetry I in the bulk.Moreover, analogously to the case of roughness-sensitive, uncompensated surface magnetization, roughness-sensitive induced surface magnetization can exist in AFMs with ðΘjtÞ preserved in the bulk, as long as there is a symmetry ðΘjt 0 Þ that combines Θ with a fractional translation t 0 perpendicular to the surface in question.
Let us now consider the form of the bulk magnetic multipoles which should correspond to this category of surface magnetization.Because we are dealing with magnetically compensated surfaces, the local-moment ME multipolization tensor must be zero; therefore, as explained in Sec.IVA, an induced surface magnetization should instead correspond to ferroically ordered atomicsite multipoles in the surface plane.The rank of the relevant multipoles, that is, whether they are atomic-site ME multipoles, magnetic octupoles, hexadecapoles, or higher rank, depends on the bulk MSG and the symmetries of the Wyckoff sites.
In contrast to the roughness-robust induced surface magnetization discussed for Cr 2 O 3 and FeF 2 in part III, however, roughness-sensitive induced surface magnetization should exist for an atomically flat surface only.This occurs if the atomic-site multipoles order ferroically within each surface plane but with opposite sign between adjacent planes.This is completely analogous to roughnesssensitive, uncompensated surface magnetization discussed in Sec.III A for which the local-moment ME multipolization, and correspondingly the surface magnetization, changes sign for multipolization corresponding to surfaces separated by atomic steps.
V. THE "NULL" CASE: ð1 10Þ NiO For (100) and ð 120Þ Cr 2 O 3 and (110) FeF 2 , we have demonstrated that, for surfaces which are magnetically compensated in the bulk AFM ground state, a symmetryallowed induced magnetization energetically stabilizes the vacuum-terminated slab.In our final example, we show with DFT calculations that inducing a finite surface magnetization is not energy lowering for a nominally compensated surface that is not FM compatible by symmetry.
We again take rocksalt AFM NiO as an example [the primitive bulk unit cell is shown in Fig. 9(a)].As a reminder, the bulk MSG is type IV C c 2=c [15.90] with preserved I as well as Θ combined with several linearly independent translations.We already showed in Sec.III B 2 that the ideal, nonreconstructed (111) surface of NiO has uncompensated but roughness-sensitive in-plane magnetization along ½11 2. We now consider the ð1 10Þ surface, with surface normal nk½1 10 that is perpendicular to both the (111) surface normal and the bulk Néel vector polarization ½11 2. Looking at an electrostatically stable slab terminated with vacuum along ½1 10 as shown in Fig. 18(a), it is clear that the ð1 10Þ surface is magnetically compensated when the surface moments retain the bulk ½11 2 Néel vector direction.A search for the operations of C c 2=c that leave the unit normal of the ð1 10Þ surface invariant modulo translations ⊥ n yields the following four operations: 1, ð2 ½1 10 j1; 1; 1Þ, ð2 0 ½1 10 j0Þ, and ðΘj1; 1; 1Þ.The corresponding surface MPG is 21 0 , which is not FM compatible because it contains Θ.Indeed, the time-reversal operation Θ reverses every component of magnetization, so there is no direction along which magnetization is symmetry allowed on the ð1 10Þ surface.
We check this by performing DFT calculations for a slab of NiO with ð1 10Þ surfaces (details are in Appendix D).We keep the Ni moments in the center two layers of the ð1 10Þ slab fixed along the bulk ½11 2 Néel vector direction, and we cant the Ni moments in the two surface layers about two orthogonal rotation axes.We first rotate about the in-plane [111] lattice vector, corresponding to an out-of-plane canting angle θ toward vacuum (θ > 0) or toward bulk (θ < 0) for both surfaces, as shown pictorially in Fig. 18(a).Separately, we cant around the out-of-plane ½1 10 surface normal, corresponding to an in-plane rotation ϕ with respect to the bulk ½11 2 direction.This is shown for the top surface with a bird's-eye view of the surface in Fig. 18(b).Both rotations induce a finite magnetization per unit area on the surface Ni layers, along ½1 10 and [111] (for θ and ϕ, respectively).
The changes in total DFT energy, given in μeV per formula unit, are shown in Figs.18(c) and 18(d) for the outof-plane canting θ and in-plane canting ϕ, respectively.In contrast to Figs. 14(a) and 14(b), where a finite surface magnetization definitively lowers the slab energy compared to the completely magnetically compensated surface, inducing a finite surface magnetization for ð1 10Þ NiO does not lower the total energy within the resolution of our DFT calculations, as expected.Note that there is an energy lowering of approximately 0.04 μeV=f:u: at θ ¼ −0.12°in Fig. 18(c), but, because it is the same order of magnitude as the penalty energy approximately 0.02 μeV=f:u:, we do not believe that it is physical.
Qualitatively, Figs.18(c) and 18(d) differ: For the in-plane canting in Fig. 18(d), the energy change is symmetric about ϕ ¼ 0, whereas the energy change for out-of-plane canting is asymmetric about θ ¼ 0 in Fig. 18(c).Specifically, the energy increase is more rapid when the canting on both surfaces is toward vacuum (θ > 0) than when the canting is toward the bulk.Different rates of energy increase for canting toward bulk versus vacuum are to be expected, since the vacuum termination on one side of the ð1 10Þ layer introduces a unidirectional anisotropy to the local environment [72].Nevertheless, it is clear that in both cases any definitive canting away from the bulk order increases the energy of the ð1 10Þ slab of NiO.We, therefore, show that the group-theory formalism to identify AFM planes with surface magnetization also correctly predicts the cases where surface magnetization is not symmetry allowed.
Finally, we note that the surfaces with no expected equilibrium surface magnetization are characterized by atomic-site ME and higher-order multipoles which are antiferroically ordered in the surface plane.Our analysis of the Ni octupoles for the ð1 10Þ surface Ni ions shown in Fig. 18 reveals that the octupolar ordering in the ð1 10Þ surface is antiferroic, consistent with the group-theory analysis.
To conclude part III of this manuscript, in Fig. 19 we give a flowchart showing the steps to categorize induced surface magnetization, given a bulk AFM and surface of interest, based on our atomic-site multipole analysis.First, one identifies the symmetry-allowed atomic-site multipoles (at all orders of the multipole expansion) in the bulk unit cell.Then, we identify the multipoles corresponding to a ME response to an electric field parallel to the surface normal.This gives us the possible directions of the induced surface magnetization.Finally, we check the ordering of the relevant multipoles on planes of atoms parallel to the surface.If all planes ferroically ordered multipoles, the induced magnetization is roughness robust.If the multipoles are ferroically ordered in plane, but their sign switches between planes, the induced surface magnetization is roughness sensitive.Lastly, if the relevant multipoles are antiferroically ordered within a single surface plane, no surface magnetization can exist even for a perfectly smooth surface.

VI. DISCUSSION AND OUTLOOK
In this work, we have provided a comprehensive classification of surface magnetization in AFMs.We have first combined and extended the previous seminal but disjoint theoretical analyses [1,2,12,15,34] to demonstrate in part I that all categories of surface magnetization, including the frequently overlooked roughness-sensitive varieties, can be described by a unifying symmetry formalism.In parts II and III, we connected this group-theory description to a universal alternative classification in terms of bulk magnetic multipoles.In particular, we showed in part II that uncompensated surface magnetization maps to local-moment ME multipoles describing bulk magnetization asymmetry at the unit-cell level, and in part III we demonstrated that more subtle, induced surface magnetization on crystallographic planes which are compensated in the bulk originate from in-plane ferroic ordering of FIG.19.Flowchart to determine, given a bulk AFM and surface of interest ðhklÞ⊥ n with no uncompensated magnetization, whether there exists roughness-sensitive or roughness-robust induced surface magnetization based on the ordering in the unit cell of the atomicsite multipoles.atomic-site multipoles describing magnetization asymmetry around individual magnetic sites.The multipole analysis for both uncompensated and induced cases also enables immediate distinction of roughness-robust and roughnesssensitive surface magnetization.step-by-step procedure to identify the existence and category of surface magnetization for any AFM surface of interest is conveniently summarized by the two flowcharts in Figs.7 and 19.Overall, we remark that in this work we exclusively focused on spin surface magnetization.At the same time, recently some progress has been made on the definition and calculation of the surface orbital magnetization [73].This is inherently more involved due to the more delicate definition of orbital magnetization within periodic boundary conditions in bulk systems [74,75].Future investigations on a possible application of the analysis presented in this work for spin surface magnetization to this new area would, thus, be of relevance.
We conclude this manuscript with some outlook.First, we emphasize that our multipole description enables the prediction of surface magnetization based on ground-state bulk properties.Bulk magnetic multipoles can be readily calculated with standard DFT codes, and several possibilities for experimentally detecting magnetic multipoles, for example, Compton and neutron scattering [21,39,76], have been recently proposed.Thus, our analysis enables the prediction of surface magnetization based on calculations or measurements of bulk properties.
A second key result is that, based on the analogy between magnetic multipoles and ME effects, all cases of surface magnetization have a symmetry correspondence to ME responses.For roughness-sensitive varieties of surface magnetization, the ME responses occurring at the isolated unit-cell level in the uncompensated case and at the atomic-site level in the induced case are antiferroic and, thus, do not result in a net ME response when considering the periodic bulk crystal.On the other hand, roughness-robust surface magnetization necessarily corresponds to a net bulk ME response for Ek n, which must be linear in E for uncompensated magnetization, and may be higher order in E for induced roughness-robust surface magnetization.In Table I, we summarize the classification of the four discussed categories of AFM surface magnetization in terms of their bulk net ME response and the corresponding type, and ordering, of bulk magnetic multipoles.
In addition to its fundamental scientific significance, our findings that roughness-robust, induced surface magnetization always corresponds to a net ME response with Ek n significantly expands the cases of surface magnetization which should be controllable by ME annealing, that is, the control of AFM domains via applied electric and magnetic fields, as mentioned briefly in the introduction.The efficacy of ME annealing for linear ME AFMs originates from a symmetry-allowed free-energy term F ∝ α ij H i E j , with H the applied magnetic field.Analogously, a net quadratic ME response M i ∝ β ijk E j E k implies a free-energy term F ∝ β ijk H i E j E k and a similar term F ∝ γ ijkl H i E j E k E l for the cubic ME response.Because every order of bulk magnetic multipole breaks time-reversal symmetry, the corresponding ME response tensor must as well, such that β, γ, and even higher-order ME response tensors all change sign for opposite signs of the bulk Néel vector.Because surface magnetization is tied to the Néel vector, this means all examples of roughness-robust surface magnetization, not just the uncompensated case corresponding to linear ME effects, should be switchable with ME annealing.
On the other hand, because roughness-sensitive surface magnetization does not correspond to a net ME effect (Table I), it cannot generally be controlled by ME annealing.However, for very thin AFM films in which the ratio of uncompensated surface magnetic dipoles to compensated dipoles in the bulk is fairly large, it should be possible to control the surface magnetization with a reasonably sized magnetic field alone.Moreover, the requirement of an atomically smooth surface to achieve macroscopic surface magnetization in the roughnesssensitive case is becoming a less formidable challenge due to improving synthesis methods, particularly in the field of thin-film deposition [77,78].Overall, we believe that further investigation of roughness-sensitive surface magnetization is highly worthwhile given the huge number of previously overlooked possibilities it opens up, in TABLE I. Classification of the four categories of AFM surface magnetization in terms of the structure of their corresponding localmoment (LM) or atomic-site (AS) multipoles and the order of the allowed net bulk ME response.RR and RS refer to roughness-robust and roughness-sensitive, respectively.For the atomic-site multipoles characterizing the induced surface magnetization, r refers to the rank of the tensor (r ¼ 2 for ME multipoles, r ¼ 3 for magnetic octupoles, etc.) which depends on the bulk AFM MSG.addition to the fact that many AFM surfaces with such magnetization are already known to exhibit exchange bias [67,79].We hope that our analysis of roughness-sensitive surface magnetization motivates its future, experimental characterization with state-of-the-art methods such as nitrogen vacancy (NV) magnetometry and spin-Hallbased magnetotransport [5,6,80].Finally, we make a few additional remarks regarding roughness-robust, induced surface magnetization separate from its connections to higher-order ME effects.Surface magnetization has recently been detected experimentally on the nominally compensated surfaces of at least three AFMs: (110) FeF 2 [13], (110) Fe 2 TeO 6 [20], and very recently in surfaces perpendicular to (001) Cr 2 O 3 [14,66].These are all consistent with our analysis that these surfaces should have roughness-robust, induced magnetization.Induced surface magnetization may also in some cases explain exchange bias at nominally compensated AFM surfaces; further investigations along these lines would be fruitful.As mentioned in Sec.IV B, surface magnetization originating from ferroic magnetic octupoles in centrosymmetric AFMs appears connected to nonrelativistic spin splitting.Overall, the recent experimental confirmations of roughness-robust, induced surface magnetization, as well as its likely ties to other intriguing spin-based phenomena, indicate that this category holds great interest from numerous perspectives.
We hope that our work motivates further studies, both experimental and theoretical, on the identification and the properties of various types of AFM surface magnetization.We also hope it spurs a comprehensive search for AFM candidates, both linear ME and otherwise, with experimentally useful surface magnetization.by assuming a semi-infinite slab with a single, vacuumterminated surface with normal n.In contrast, it is apparent, for example, from Figs. 14(a) and 14(b) for Cr 2 O 3 the electrostatically stable, symmetric ð 120Þ and (100) slabs have two surfaces and, hence, actually belong to higher-symmetry point groups than the MPGs we analyze, which correspond to semi-infinite slabs.This is also the case for the DFT slabs we use for (110) FeF 2 and ð1 10Þ NiO.In fact, the point groups of the symmetric slabs are not even FM compatible, in general.However, we argue that using such a slab to explore surface magnetization is not problematic as long as the number of bulklike layers between the two surfaces is sufficiently large: Then, the surfaces behave independently as though they were on two oppositely oriented, semi-infinite slabs.From an energetic perspective, only the local symmetry, i.e., the symmetry of one vacuum-terminated surface and the nearby bulk with which it interacts, is relevant in determining whether an induced magnetization is favorable for this surface.And this local surface symmetry is identical to the surface MPG determined by the grouptheory and multipole formalisms.
To confirm that the local symmetry, rather than the symmetry of the entire DFT slab, is the relevant MPG to consider, we take the example of ð 120Þ Cr2O 3 and force the slab with two surfaces to have the same MPG as the idealized semi-infinite slab with a single [010]-oriented surface normal.We do this by substituting the bottom Cr layer with a monolayer of Fe as shown in the top in Fig. 20, thereby breaking the extra mirror plane and inversion symmetries created by having two equivalent surfaces in Fig. 14(a) in the main text.Because Fe naturally adopts a þ3 valence state, this "semi-infinite" slab is electrostatically stable and there are no convergence issues.For this calculation, we use the same Hubbard U of 4 eV for both Fe and Cr 3d states.We use the same convergence criteria as for the symmetric ð 120Þ slab and use the same Lagrange multiplier λ ¼ 10 for the constrained magnetic calculations.
We fix the Fe moments as well as the Cr moments in the center two layers along [001] in the "out" domain, and we again plot the change in total energy as a function of canting angle for the top layer of Cr.While there are small quantitative differences between Figs. 20 and 14(a), the results for the ð 120Þ surface with an Fe monolayer in Fig. 20 exhibit the same energy lowering with a minimum at a 0.25°canting angle toward the vacuum.This demonstrates that, indeed, it is the local surface symmetry for the DFT slab, reflecting the semi-infinite slab MPG, which dictates the energetics of the induced surface magnetization.Therefore, for all DFT calculations in the main text, we use the electrostatically stable symmetric slabs, which have higher global symmetry than the local surface MPG, to reduce computational cost.For density-functional calculations for the (110) FeF 2 surface, we again use VASP.We employ the generalized gradient approximation using the Perdew-Burke-Ernzerhof (PBE) functional [87], and we select a Hubbard U of 6 eV and a Hund's exchange J of 0.95 eV.These parameters are based on previous DFT þ U calculations of FeF 2 and its (110) surface [68,69].Because surface magnetization in the case of FeF 2 develops collinearly along the bulk Néel vector in contrast to ð 120Þ Cr 2 O 3 , we do not include SOC in our calculations and instead use collinear spin-polarized DFT.In calculating the site-projected magnetization for the Fe and F ions, we use the default integration spheres of 1.058 and 0.794 Å, respectively.We relax the bulk FeF 2 structure using the plane-wave and force criteria described for Cr 2 O 3 in the previous section and a 6 × 6 × 9 Gammacentered k-point mesh for the bulk tetragonal cell.We then create a electrostatically stable (110) slab with eight Fe layers and 15 Å of vacuum along ½110.Note that we increase the slab thickness of (110) FeF 2 compared to ð 120Þ Cr 2 O 3 due to previous reports that approximately eight layers are required to obtain bulklike behavior in the center of the slab [68,69], which was not the case for ð 120Þ Cr 2 O 3 .Finally, we relax all internal coordinates of FIG.20.Top: slab of ð 120Þ chromia with a monolayer of Fe replacing the bottom Cr, such that the MPG of the total slab is lowered and matches that of a semi-infinite slab (MPG 2).Bottom: total DFT energy per formula unit for the Cr-terminated surface as a function of the canting angle.The penalty energy E p in the constrained magnetic calculations is plotted with a dashed line and star markers.
the slab while keeping the lattice vectors fixed to our bulk relaxed values.
For the calculation of the lattice-mediated ME response, phonon frequencies and eigendisplacements are computed using density-functional perturbation theory [88], as implemented in the QUANTUM (QE) [89,90] and THERMO_PW [91] packages.Self-consistent ground-state and linear response calculations are performed within the collinear spin-polarized generalized gradient approximation scheme, without SOC included, and the exchangecorrelation energy is treated using the PBE parametrization [87].Ion cores are described using scalar-relativistic ultrasoft pseudopotentials (PPs) [92], with 4s and 3d valence electrons for Fe (PP Fe. pbe-n-rrkjus_psl.1.0.0.UPF from pslibrary 1.0.0 [93,94]) and with 2s and 2p valence electrons for F (PP F. pbe-n-rrkjus_psl.1.0.0.UPF).Correlation effects are dealt with by introducing a Hubbard U correction [95], with U ¼ 6 eV and J ¼ 0.95 eV, selected to be consistent with VASP calculations on the (110) surface.Pseudo wave functions and charge density are expanded in a plane-wave basis set with kinetic energy cutoffs of 80 and 400 Ry, respectively.BZ integrations are performed using a Γcentered uniform Monkhorst-Pack [85] mesh of k points with 4 × 4 × 6 points.
Frozen phonon calculations of the induced Fe magnetic moments are performed at the distorted structure obtained by freezing in the phonon ionic displacements.The induced local magnetic moments are computed using the allelectron linearized augmented-plane-wave (LAPW) code ELK [96], to have a more accurate description of the atomic magnetic moments, within the local spin-density approximation, with spin-orbit coupling included.The APW wave functions and potential inside the muffin-tin spheres are expanded in a spherical harmonic basis set, with cutoff l max ¼ 12. Similarly to the calculations performed in QE, a Hubbard U correction with the same values of U and J is applied.The BZ is sampled with the same mesh of k points used in QE.

APPENDIX D: DFT CALCULATION DETAILS FOR ð1 10Þ NiO
For density-functional calculations of ð1 10Þ NiO, we again use VASP and choose the generalized gradient approximation with the PBE functional [87].We take U ¼ 4.6 eV within the Dudarev implementation [84] of DFT þ U [83].This is based on the parameters used for DFT þ U calculations of NiO in Ref. [49].We fully relax the bulk structure, and then relax the internal coordinates of a four-layer slab with an electrostatically stable ð1 10Þ surface and 15 Å of vacuum along ½1 10.We note that the slab is quite thin (approximately 6 Å); due to the large number of atoms in the ð 110Þ surface plane [96 for the slab shown in Figs.18(a) and 18(b)], a thicker slab would be computationally prohibitive.By checking the layerprojected density of states, we confirm that the center of the slab is sufficiently bulklike.Analogously to the case of ð 120Þ Cr 2 O 3 , we perform constrained magnetic calculations with a Lagrange multiplier λ ¼ 10.

APPENDIX E: MAGNETIC OCTUPOLE TENSOR
Here, we report the spherical irreducible components of the magnetic octupole tensor O ijk introduced in Sec.IVA.For the full derivation, see Ref. [22].O can be decomposed into a totally symmetric tensor S and a residue tensor R: The traceless part T of the totally symmetric tensor S is expressed in terms of seven independent parameters (O −3 , O −2 , …, O 2 , O 3 ) as follows: ; We refer to these parameters as proper magnetic octupoles, because they transform into each other upon rotations in the same way as the atomic f orbitals.
The residual tensor R can be further decomposed into two tensors, written by using five and three independent parameters, respectively.The five-dimensional residue tensor R ð5Þ is written in terms of the quadrupole moments Q ðτÞ of the toroidal moment density τðrÞ ¼ r × μðrÞ as

FIG. 3 .
FIG. 3. "Roughness-robust" versus "roughness-sensitive" surface magnetization.(a) Roughness-robust: any operation in the bulk MSG preserving the direction of n and containing a translation component k n leaves m FM ðg nÞ invariant.(b) Roughness-sensitive: at least one operation in G with t= ⊥ n reverses the direction of m FM ðg nÞ.

FIG. 6 .
FIG. 6. Uncompensated surface magnetization for two different terminations of (001) Cr 2 O 3 calculated from the local-moment multipolization.(a) Left: bulk unit cell defining the electrostatically stable termination of (001) chromia (right).The [001] oriented surface magnetization for this termination is 11.9μ B =nm 2 from the ðz; zÞ component of MLM .The polar, oxygen-terminated unit cell shown in (b) is related to this unit cell by translating the topmost Cr atom down by one lattice vector, shown by the blue arrow and faded-out atom at the bottom of the cell.(b) Left: bulk unit cell related to the cell in (a) by one multipolization increment.This basis defines the polar (001) surface (right) with surface magnetization −2.4μ B =nm 2 .

FIG. 8 .
FIG. 8. (a) Primitive rhombohedral cell and magnetic ordering of Cr 2 O 3 .Cr and O atoms are identified by blue and red spheres, respectively.(b) Top view of the conventional hexagonal cell of Cr 2 O 3 .
(a).The ground-state AFM order consists of FM planes stacked antiferromagnetically along the [111] direction, with the Ni magnetic moments polarized in the planes along the ½11 2 and ½ 1 1 2 directions [49].NiO has the bulk type-IV MSG C c 2=c [15.90].When discussing the symmetries of NiO, it is useful to adopt the basis of a magnetic supercell in which the Ni moments are perpendicular to two of the lattice vectors and parallel to the third.To this end, we use a magnetic cell with lattice vectors a ¼ ½11 2, b ¼ ½1 10, and c ¼ ½222 in terms of the lattice vectors of the conventional nonmagnetic cell in Fig. 9(a).The resulting, reoriented supercell is shown in Fig. 9(b) for two opposite AFM domains.Now let us consider the (111) surface of NiO, which is perpendicular to the direction of AFM stacking.The rotated unit cells shown in Fig. 9(b) show an atomically smooth termination of the (111) surface, with the two AFM domains having opposite orientations of the Ni moments in the topmost layer.Note that, due to the alternating, equally spaced layers of Ni 2þ and O 2− ions along [111], it is not possible to create an atomically smooth electrostatically stable (111) surface.Both Ni and O terminations are intrinsically polar, and, thus, the true (111) surface of NiO almost certainly undergoes a reconstruction [58,59].We nevertheless consider the unreconstructed, polar (111) surface here as a straightforward example of a Θ-symmetric AFM surface with uncompensated roughness-sensitive magnetization.Inspection of either structure in Fig. 9(b) shows that such a surface has an uncompensated FM layer of Ni moments along either ½11 2 or ½ 1 1 2.

FIG. 9 .
FIG. 9. Uncompensated, roughness-sensitive surface magnetization in (111) NiO and (001) Fe 2 O 3 .(a) Nonmagnetic unit cell of rocksalt NiO.(b) Unreconstructed (111) surface of NiO.Left and right surfaces, with equal magnitudes and opposite signs of surface magnetization based on the local-moment multipolization, are connected by time reversal plus a translation along the [111] surface normal (dashed black arrow).(111) NiO has a "half-increment-containing" multipolization array.(c) Two inequivalent atomic terminations for (001) Fe 2 O 3 .The left two figures correspond to a polar O termination, and symmetrically equivalent surfaces with equal magnitudes and opposite signs of surface magnetization are related by a mirror plane containing [001], plus a fractional translation along [001].Right: Two symmetrically equivalent surfaces with the electrostatically stable termination, which has zero surface magnetization.(001) Fe 2 O 3 has a "zero-containing" multipolization array.
)] about [001] as for the R 30 c 0 MSG of Cr 2 O 3 , so the surface MPG of (001) Fe 2 O 3 , 3, is FM compatible and can host surface magnetization along [001].

2 . 1 .
Induced surface magnetization in ð 120Þ and (100) Cr 2 O 3 We first look again at ME AFM Cr 2 O 3 , with broken I and Θ in the bulk MSG R 30 c 0 .Bulk Cr 2 O 3 can be described by a primitive rhombohedral cell [see Fig. 8(a), where we show also the magnetic dipolar order of Cr atoms] or likewise by a conventional hexagonal cell, whose top view is shown in Fig. 8(b).In contrast to the uncompensated (001) surface we analyze in Sec.III B 1, this time we examine a surface perpendicular to the (001) plane.In Fig. 10(a), we indicate the high-symmetry directions in the hexagonal basal plane.Note that, for surfaces perpendicular to (001) in hexagonal Cr 2 O 3 , the real-space surface normals and their corresponding Miller planes do not have the same indices.ð 120Þ surface, symmetry analysis.We focus first on the ð 120Þ surface.The bulk, orthorhombic unit cell which defines the electrostatically stable ð 120Þ surface when tiled semi-infinitely along the perpendicular [010] direction is shown in Fig. 11(a).The magnetic lattice consists of two layers with 12 Cr each stacked along [010].Note that, within each layer, half of the Cr magnetic moments are oriented along [001] in the ground state, and half are pointed along ½00 1.Therefore, any flat ð 120Þ surface will be magnetically compensated.
FIG. 10.(a) High-symmetry directions in the basal plane of a hexagonal cell.(b) Symmetry operations of the 3m point group.

FIG. 11 .
FIG. 11.Bulk unit cells which can be tiled semi-infinitely to define the electrostatically stable atomic terminations of the (a) ð 120Þ and (b) (100) surfaces of Cr 2 O 3 .The magnetic dipoles moments of the Cr ions, shown in turquoise, are ordered according to the bulk AFM ground state with the Néel vector parallel to [001].
FIG.12.Third-order ME response γ zxxx to an electric field pointing along any direction in the basal plane of the conventional hexagonal cell.θ identifies the angle between the direction of the electric field and the x axis.

FIG. 13
FIG. 13.(a) Top view of the ð 120Þ Cr 2 O 3 surface for the outpointing (top) and in-pointing (bottom) domain.(b) Side view of the (100) Cr 2 O 3 surface for the out-pointing (top) and in-pointing (bottom) domain.

FIG. 14 .
FIG. 14.Energy ΔE per formula unit as a function of the canting angle for the in-pointing (green) and out-pointing (red) domain, and net out-of-plane surface magnetization as a function of the canting angle for (a) ð 120Þ and (b) (100) Cr 2 O 3 surfaces.Penalty energies for in and out domains are plotted with green circles and red stars, respectively.Dashed lines in the bottom of both (a) and (b) identify the surface magnetization corresponding to the canting angle that minimizes the energy.Note that, for the (100) surface, an in-plane component of surface magnetization along [001] is also allowed to develop; we depict the difference in magnitudes of the [001] components of sublattice magnetization by the different colored arrows for the surface moments in the cartoon of the slab.

FIG. 15 .FIG. 16 .
FIG. 15.(a) Unit cell and magnetic ordering of FeF 2 .Fe and F atoms are identified by gold and gray spheres, respectively.(b) Top view of the unit cell of FeF 2 , with FeF 6 octahedra highlighted.

FIG. 18 .
FIG. 18. Induced magnetization on the non-FM-compatible ð1 10Þ surface of NiO.(a) Electrostatically stable slab used in our DFT þ U calculations for the ð1 10Þ surface, demonstrating the out-of-plane rotation θ for the surface Ni moments (in orange) which we enforce in our constrained magnetic calculations.Atoms deeper into the page are removed for ease of visualization.(b) View of the slab looking down at the ð1 10Þ surface, showing the in-plane rotation ϕ for the surface moments.Faded-out blue moments are in the layer below the surface.(c) and (d) show the changes compared to the energy with the bulk ½11 2 polarization as a function of the rotation angles in (a) and (b), respectively.Blue lines in (c) and (d) show the penalty energy E p in the constrained magnetic calculations.
Ferroic M AS r in-plane, antiferroic along n Antiferroic M AS r in-plane Bulk ME responses for Ek n Linear M ∝ E None M ∝ E r−1 None None APPENDIX C: DFT CALCULATION DETAILS FOR (110) FeF 2 Tables summarizing the categories of surface magnetization that can exist for certain Miller planes in a given bulk AFM.In the left table, we again depict cartoons for the four categories of surface magnetization we discuss.In the right table, we list the possible categories of surface magnetization based on the absence or presence of Θ and I in the bulk MSG of the AFM of interest.ðΘjtÞ represents time-reversal Θ combined with any fractional translation t.