{'fraction_non_alphanumeric': 0.03227221026955175, 'fraction_numerical': 0.007153330683951599, 'mean_word_length': 4.359036498494529, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Towards a theory of surface orbital magnetization

The theory of bulk orbital magnetization has been formulated both in reciprocal space based on Berry curvature and related quantities, and in real space in terms of the spatial average of a local marker. Here we consider a three-dimensional antiferromagnetic material having a vanishing bulk but a nonzero surface orbital magnetization. We ask whether the surface-normal component of the surface magnetization is well defined, and if so, how to compute it. As the physical observable corresponding to this quantity, we identify the macroscopic current running along a hinge shared by two facets. However, the hinge current only constrains the difference of the surface magnetizations on the adjoined facets, leaving a potential ambiguity. By performing a symmetry analysis, we find that only crystals exhibiting a pseudoscalar symmetry admit well-defined magnetizations at their surfaces. We then explore the possibility of computing surface magnetization via a coarse-graining procedure applied to a local marker. We show that multiple expressions for the local marker exist, and apply constraints to filter out potentially meaningful candidates. Using several tight-binding models as our theoretical test bed and several potential markers, we compute surface magnetizations for slab geometries and compare their predictions with explicit calculations of the macroscopic hinge currents of rod geometries. We find that only a particular form of the marker consistently predicts the correct hinge currents.


I. INTRODUCTION
The modern theories of bulk electric polarization P and bulk orbital magnetization M [1][2][3][4][5][6][7] express these quantities as Brillouin zone integrals of quantities involving Berry connections and curvatures of ground state Bloch functions. However, the ground-state of systems of independent electrons may also be uniquely described by the single-particle density matrix, also known as the ground state projector P (r, r ). This is a quantity that for insulating states of matter decays exponentially with |r − r |, even for Chern insulators [8].
A natural question to ask is whether it is possible to express P and M in terms of P (r, r ). In the case of the polarization, while the change ∆P can be determined by following the change of P (r, r ) between two states connected by an adiabatic switching process, there is no corresponding expression for P itself, which is determined only modulo a quantum [1,2].
However, M does not suffer from any quantum of indeterminacy, and should be expressible via the ground state projector. Bianco and Resta [9] demonstrated that this is, in fact, correct. Starting from an expression for a simple 2D crystallite of finite area, the authors demonstrated that for any insulator, even a Chern insulator, the orbital magnetization M is given in the thermodynamic limit by Here A is the bulk unit cell area and M(r) is a local function, hereafter referred to as the local marker, which is expressed in terms of P (r, r ). Although M(r) appears to play the role of magnetic dipole density, only it's macroscopic average bears any physical meaning. This is in keeping with the fact that while microscopic charge and current densities are well defined, microscopic dipole densities (be it charge or orbital magnetic) are not [10]. Recent years have witnessed numerous efforts aimed at defining and providing explicit expressions for surfacespecific analogues of bulk quantities. Examples include the surface anomalous Hall conductivity (AHC) [11,12] and the surface-parallel boundary electric polarization [13][14][15][16][17] on the surface of a bulk insulator. To date, however, there has been very little discussion of surface orbital magnetization in the literature [18], and a proper theory of orbital magnetization at the surface of a bulk material has yet to be developed.
In the present work, we explore the possibility of defining surface orbital magnetization. By this we mean the excess surface-normal macroscopic magnetization (magnetic moment per unit area) at the surface of a bulk material with broken time-reversal (TR) symmetry. Our investigation is restricted to antiferromagnetic systems featuring a vanishing bulk magnetization, as this allows us to readily disentangle surface contributions to the orbital magnetization from those stemming from the bulk. We shall also restrict ourselves throughout to the case of insulating surfaces of insulating crystals, leaving aside for now any special considerations that might arise for metallic surfaces.
The case of surface spin magnetization is much more straightforward, as the latter can be obtained simply by integrating the net spin density ∆n = n ↑ − n ↓ over the surface region after an appropriate coarse-graining (e.g., using window-averaging methods). In the remainder of this work, therefore, we focus solely on the orbital component of the surface magnetization. Here it is natural to expect more difficulties, since even the theory of bulk orbital magnetization has been put on a solid footing only in relatively recent times [3][4][5][6][7]. The essential problem, as in the theory of electric polarization, is that the position operator is ill defined in the Bloch representation, so that methods based on Berry connections and curvatures are required instead. Another way to frame the problem is to note that expressions are available for neither the local polarization density P(r) nor the local orbital magnetization density M(r).
Even from a phenomenological point of view, another difficulty arises. The physical observable corresponding to the presence of surface magnetization is the presence of macroscopic bound currents running along the hinges where surface facets meet, with the hinge current given by the difference between the surface magnetizations on the two adjoining facets. However, since the hinge currents only determine surface magnetization differences, this raises the question whether surface orbital magnetization can be uniquely defined, or only defined up to a constant shift in the values predicted for all facets of a given bulk material. Such a "shift freedom" is not easily ruled out a priori.
In this work, we first investigate the symmetry conditions that allow such a shift freedom to be ruled out, so that surface magnetizations are clearly uniquely defined. We then explore whether the use of a local marker M(r), such as the one proposed by Bianco and Resta, can be adopted for the purpose of defining a surface orbital magnetization. That is, even if the value of M(r) at a single point has no obvious physical meaning, we ask whether it is possible to coarse-grain and integrate such a marker over the surface region to obtain a valid expression for the surface orbital magnetization. This paper is organized as follows. In Sec. II, we identify the physical observable corresponding to the presence of a surface magnetization, namely the macroscopic current running along a hinge shared by adjacent surface facets of a bulk crystal, given by the difference between the magnetizations of the two surface facets forming the hinge. We show that this single relation is insufficient to ascribe an unambiguous value of magnetization to any of the facets; constraints in the form of crystalline symmetries are required to accomplish this. In Sec. III, we identify all the possible symmetries that lead to unambiguous surface magnetizations. We find that this set of symmetries is identical to the set of symmetries that quantizes the Chern-Simons axion coupling.
In Sec. IV we introduce in greater detail the local marker formulation of orbital magnetism, and introduce our formalism for the calculation of surface magnetization and hinge currents. The results of Sec. III and Sec. IV are later tested in Sec. V by calculations performed on tight-binding models. We end our paper with a discussion of our results and some speculations on their connection to the theory of orbital magnetic quadrupole moments in Sec. VI, and summarize in Sec. VII.

II. SURFACE MAGNETIZATIONS AND HINGE CURRENTS
Consider a stand-alone 2D system such as a monolayer or multilayer with layer-normal magnetization M ⊥ . The magnetization manifests itself on any edge of the system as a bound current of magnitude M ⊥ , with its sign determined from the right-hand rule. If instead we consider the surface of a 3D bulk system, we would also expect any surface-normal magnetization at the surface to manifest itself at an edge in a similar fashion. However, in the presence of the bulk, the objects playing the role of edges are the hinges, i.e., the intersections of neighboring surface facets. So, while in the case of an isolated 2D system the edge current directly determines the magnetization M ⊥ , any current present on a hinge only determines the difference of M ⊥ values on the two surface facets meeting at the hinge. We depict this situation schematically in Fig. 1, where it is evident that the hinge current I LR is given by Note that a line current may also be present at a line defect, such as a step or domain wall, that separates surface patches having the same orientation but different M ⊥ values. We can regard these as 180 • hinges, and the treatment of such cases introduces no new difficulties. In the remainder of this work, we restrict the discussion to true hinges connecting facets of different orientation. In any case, we emphasize that the line currents in question are bound currents localized at insulating hinges, not the chiral conductance channels that may occur in topological systems, e.g., at the boundary of a 2D quantum Hall state or on the hinges of an axion insulator [19,20].
For the entirety of this paper, we take the hinge current to be a physical observable, being detectable in principle by the circulating magnetic field it creates. Combined with Eq. (2), this assumption indicates that differences between magnetizations on neighboring surface facets are observables. However, there is no a priori reason to believe that the individual values of facet magnetizations are well defined. We must therefore understand whether, or under what circumstances, we can resolve this ambiguity and assign a unique M ⊥ to a given facet of the bulk. Once an unambiguous value of magnetization is assigned to even a single facet, we may then repeatedly employ Eq. (2) to find the magnetization at any other facet of the 3D crystallite.
While we have focused above on the surface-normal component of the surface orbital magnetization, an inplane component M could also be present. If so, its strength is easily computed from a detailed knowledge of the microscopic current distribution in the surface region, unlike the more problematic M ⊥ . In principle, M should where S and N refer ton-symmorphic andn-nonsymmorphic operations respectively. Superscripts '−' and '+' indicate whether or not the operation reverses the sign of the outwarddirected M ⊥ on a surface with normaln, which also corresponds to whether the symmetry operation is axion-odd or axion-even.
be observable by the presence of a corresponding concentration of B in the surface region (since H vanishes in this geometry in the absence of free current). When M is present, we can filter out its contribution to the hinge current by constructing an appropriate Amperian loop that is normal to both surface facets, as depicted in Fig. 2. Integrating the magnetization around this loop yields the bound current passing through the loop. Since we assume the absence of free current, and since the vacuum and bulk are free of magnetization, only the top and side M ⊥ values contribute to the hinge current. We adopt this as the definition of the hinge current in such situations.

III. SYMMETRIES AND SURFACE MAGNETIZATION
The goal of this section is to establish under what circumstances we can define an unambiguous magnetization on a given surface facet of a 3D crystallite. As noted in the previous section, Eq. (2) effectively results in one equation in two unknowns, which prevents an assignment of definite values to M L ⊥ and M R ⊥ . We therefore need additional restrictions on the possible values M L ⊥ and M R ⊥ may take. These restrictions will appear in the form of symmetry-induced constraints. We first introduce the notation we will be using to describe the actions of the symmetry operations, and subsequently delve into the symmetry analysis of surface magnetization.

A. Notation
Consider a surface with unit normaln. Let the outward surface-normal component of the surface magnetization be M ⊥ = M·n, with units of magnetic moment per unit area. With respect to this choice ofn, we will focus on operations that either preserve or reverse the direction n, as well as those that preserve or flip the sign of M ⊥ . Operations that flipn will be denoted by R ± n ; the superscript '+' or '−' indicates that M ⊥ is preserved or flipped, respectively. We further classify operations that preservê n as either symmorphic or nonsymmorphic along that direction. Thus,n-symmorphic operations will be denoted S ± n (even if they involve fractional translations in the plane normal ton), whilen-nonsymmorphic operations will be denoted N ± n . We summarize the classification of symmetries in Table  I. We note that operations labeled with a superscript '−' are axion-odd while those labeled with a '+' superscript are axion-even. By definition, an axion-odd symmetry is either a proper spatial rotation with time reversal or an improper rotation without time reversal. The significance of this distinction will become clear in Sec. III C.
For concreteness, letn =ẑ. We list all the symmetry operations that fall into the categories listed above in Table II. In our notation, E denotes the identity operation and I denotes inversion. C n denotes a rotation by 2π/n about theẑ axis (n = 2, 3, 4, 6) while C 2 is a 2-fold rotation about an axis in the x-y plane. m z denotes a mirror reflection across the x-y plane, and m d is a mirror reflection across a plane containing theẑ axis. S n = m z C n denotes an improper rotation (rotoinversion) about thê z axis. A prime indicates composition with time reversal. Nonsymmorphic fractional translations alongẑ are indicated as c/2 or pc/n, where c is the vertical lattice constant and p < n is an integer (or half-integer in the case of C 3 ). Half lattice translations in the x-y plane are indicated by τ d/2 . In what follows, we shall often denote the screw operation {C n |pc/n} as just n p for conciseness. It is important to note that some of the above operations may belong to R, S, and N operations with respect to other directions. For example, the glide mirror {m y |c/2} is both an N − z and an R − y operation.

B. Symmetry Analysis for Surface Magnetization
In this subsection, we explore how to exploit the symmetry operations in the S ± n , R ± n , and N ± n classes to define an unambiguous surface magnetization on at least one surface facet; repeated application of Eq. (2) will subsequently help define magnetizations for all other surface facets. The main idea is that a given symmetry operation will relate specific surface facets to one another, allowing us to explicitly relate physical quantities defined on those facets. We then consider the magnetizations M ⊥ at the surfaces, and check whether this configuration allows for an unambiguous determination of these M ⊥ values. We find that S − n , R − n , and N − n operations allow us to define unambiguous surface magnetizations, while for S + n , R + n , and N + n this is not possible. As a reminder, we are focusing onn =ẑ for concreteness. We first investigate the axion-odd symmetries S − z , R − z , and N − z , and subsequently the axion-even symmetries S + z , R + z , and N + z .
FIG. 3. Schematic of a slab geometry for a bulk material with an R − z symmetry. An additional facet of arbitrary orientation appears at right. The R − z symmetry ensures that the magnetizations on the top and bottom facets are equal and parallel, as indicated by the maroon arrows. All three facet magnetizations may be found using the known hinge currents I1 and I2, as discussed in the text.
Symmetries of type S − z . Suppose we find a symmetry of the type S − z in the bulk symmetry group. We may then construct a surface respecting this symmetry with unit normalẑ. Under an S − z operation, the surface will remain invariant, but the magnetization at the surface will be flipped. This can only mean that the magnetization at the surface is zero.
Symmetries of type R − z . If we find an R − z symmetry in the bulk symmetry group, we may construct a finite slab in z consistent with the identified R − z symmetry applied to the slab as a whole. This construction will then enforce identical magnetizations on top and bottom surfaces.
Let us introduce an arbitrary third surface facet to the slab, as shown in Fig. 3. In the figure, the arrows normal to the surfaces indicate the direction of the surface-normal magnetization in the given global Cartesian frame. The vectors are labeled by their corresponding magnitudes; since M B ⊥ < 0, the magnitude of the bottom surface magnetization is labeled as −M B ⊥ . We then also have two hinges with currents I 1 and I 2 , with the directions of the arrowheads on the hinges indicating the directions of positive current flow.
Using our knowledge of these currents and the symmetry of the bulk slab, we may then write down a system of equations  {m d |c/2} and {C 2z |c/2} are also R − d operations. We therefore already know that we can define unambiguous surface magnetizations in those cases. We thus focus on the time-reversed screws with 3-fold, 4-fold, and 6-fold rotation axes.
For these operations, we construct crystallites in the shape of regular triangular, square, or hexagonal prisms, respectively, as shown schematically in Fig. 4. Each prism is infinite along z and is constructed to respect the corresponding screw symmetry applied to the rod as a whole. In the case of 3 m , symmetry dictates that each facet have vanishing magnetization. For the 4 m and 6 m screws, the surface magnetizations may be found from the currents on the prism hinges.
Symmetries of type S + z . We construct a surface as in the S − z case. Although S + z operations leave the surface invariant, they do not flip M ⊥ . We therefore cannot extract any meaningful information from this construction.
However, operations in this class may also relate surfaces with normals not alongẑ. For the case of C n operations with n = 3, 4, 6, we may generate the same infinite prism constructions as in the N − z class of operations. However, as we see in Fig. 5, all the hinge currents vanish, and although we know the relative orientations of the facet magnetizations, we cannot extract any information on their magnitudes as a result.
For the C 2 as well as m d operations, we may construct a slab geometry as in the case of R − z operations, except that the z-direction will be periodic, and the slab is consistent with either the C 2 or the m d symmetry. As a result, any possible surface magnetizations will be equal and antiparallel on the two surfaces, as shown for the particular case of an m x symmetry in Fig. 6. With the introduction of an additional arbitrary surface facet, we may attempt to solve a system of equations Symmetry constraints arising from 3m, 4m, and 6m symmetries. Conventions are the same as in Fig. 4. (a) Triangular prism constructed using 3m symmetry. (b) Square prism constructed using 4m symmetry. (c) Hexagonal prism constructed using 6m symmetry. In each case all facets have the same surface magnetization M ⊥ , whose magnitude cannot be determined since hinge currents are absent.
involving the hinge currents. However, this system is unsolvable. Thus, S + z operations alone are unable to unambiguously define surface magnetizations.
Symmetries of type N + z . Being nonsymmorphic operations in theẑ direction, no surface with the same unit normal will ever be preserved by an N + z operation, so these operations by themselves tell us nothing about magnetizations at this surface. Furthermore, they provide no information on surfaces transverse to theẑ surface, as the arguments in Figs. 5 and 6 apply.
Symmetries of type R + z . The same arguments as presented for the S + z and N + z cases hold here as well. We therefore cannot define surface magnetizations in this case either.

C. Discussion
Our symmetry analysis reveals that we may unambiguously define magnetizations on the surfaces only when the bulk symmetry group features S − n , R − n , or N − n operations.
In these cases, the argument proceeded via a thought experiment involving the construction of a slab or rod in such a way that the surface facets are related to each other by a global application of the relevant bulk symmetry. We emphasize that, for these bulk symmetries, the surface magnetizations are still well defined even for configurations in which different facets are terminated in different ways. The conclusion that M ⊥ is well defined follows just from knowing that it is possible, in principle, to prepare a globally symmetric slab or rod.
As it happens, the S − n , R − n , and N − n symmetries that allow an unambiguous definition are precisely the ones that are classified as axion-odd according to Table I. In other words, these are the symmetries that reverse the sign of the Chern-Simons axion (CSA) coupling, which describes a topological contribution to the isotropic linear magnetoelectric coupling of the material [12,21,22]. This is characterized by an "axion angle" θ CS that takes the form where the integration is performed over the Brillouin zone and A mn µ (k) = i u mk |∂ µ u nk is the Berry connection matrix, with |u nk being the cell-periodic part of the Bloch function of occupied band n.
A defining feature of the CSA coupling is that it is gauge-invariant modulo 2π under a unitary mixing of the occupied ground state Bloch functions. That is, if we perform the transformation where U mn (k) is a unitary matrix, then θ CS → θ CS +2πn where n is an integer. The operations of the type S − n , R − n , and N − n reverse the sign of θ CS ; but since the coupling is defined modulo 2π, it follows that θ CS = 0 or π. Thus, the axion angle is quantized to values of 0 or π in the presence of axion-odd symmetries. By contrast, an axion-even symmetry (simple rotation or time-reversed improper rotation) allows for a generically nonzero value of θ CS .
This identification of an intimate connection between the ability to define unambiguous surface magnetizations and the presence of a quantized axion coupling is one of the principle contributions of this work.
There is another suggestive connection between the physics of the axion coupling and the surface magnetization. The orbital magnetization of a 2D insulator is known to have a dependence on the Fermi energy µ, even for µ in the gap [5,9,23], given by the Středa relation dM/dµ = σ AHC /e, if an anomalous Hall conductivity σ AHC is present. We also know that the surface σ AHC of an insulating surface of a 3D insulator is constrained, modulo e 2 /h, by the axion coupling [11,12,20,23]. The application of these principles to our case suggests that, for a crystal having a generic nonzero axion coupling and a common branch choice of the quantum on all surfaces, 1 the M ⊥ values of all facets will undergo a common constant shift in response to a global shift of the Fermi level in the gap. This is already a kind of a realization of the "shift freedom" discussed in Sec. I.
While we have been focusing on the consequences of symmetry in this section, this can only take us so far. Even in the axion-odd case where we expect a unique definition of surface magnetization, the symmetry arguments do not answer the question of how, in principle, this quantity can be obtained from a direct calculation using only information about the surface Hamiltonian and surface electronic structure of a given surface. In the opposite case that no axion-odd symmetry is present, the same question arises, but now there is the additional issue of a possible shift freedom in the definition of M ⊥ . To go beyond the symmetry arguments, we need to develop an explicit framework for carrying out calculations of surface orbital magnetization in the general case. The remainder of this work is devoted to an exploration of the use of orbital magnetization local markers for the construction of such a framework.

IV. METHODS
We now turn to a discussion of the computational aspects of surface magnetization and hinge currents. We expand on the local marker formulation of orbital magnetization [9] and discuss its application in computing surface magnetization. The section is concluded with a discussion of microscopic currents and an explanation of their relationship to the macroscopic hinge current.

A. Local-marker formulation
We now work in 3D, where Eq. (1) takes the form Here V is the unit cell volume and both M and M(r) have units of magnetic moment per unit volume. (Quantities denoted as M ⊥ or M surf will continue to be 2D surface magnetizations with units of magnetic moment per unit area.) We focus initially on the single component M z of M, dropping the z subscript for conciseness. The local-marker formulation is a direct result of expressing the orbital magnetization in terms of the single particle density matrix P (r, r ). Bianco and Resta [9] demonstrated that for bulk insulating systems at T = 0 at the single-particle level, where P 2 = P , the magnetization M may be written as with where e < 0 is the electron charge, 2 the trace is taken per unit cell of volume V , Q = 1 − P , H is the Hamiltonian, and µ is the Fermi energy. M LC , M IC , and M C correspond to the so-called local circulation, itinerant circulation, and Chern number contributions to the magnetization, respectively [5]. Eq. (7), expressing the total M as a trace of a local marker M(r) over position space, is evidently satisfied by choosing with M IC (r) = e Im r|QxP HP yQ|r , Note that M C (r) is directly proportional to the local Chern marker C(r) [11,24], i.e., consistent with the discussion towards the end of Sec. III C.
It is important to note that due to the invariance of the trace under cyclic permutations of operators, different expressions for M(r) stemming from Eq. (8) are possible. For example, M LC (r) in Eq. (12) could be replaced by achieved by applying a cyclic permutation to M LC (r) in Eq. (13) and using the fact that Q 2 = Q. The L and R superscripts indicate that H has been moved to the "left" or "right" respectively, and henceforth we shall refer to the marker in Eq. (13) as M C LC (C for "center"). Integrating Eq. (17) or (18) over all space yields Eq. (9), just as integrating Eq. (13) does. We refer to the different expressions for the local marker as trace projections. Such an ambiguity in the local marker may be viewed as a manifestation of the fact that microscopic magnetic dipole densities are not well defined.
We can argue that there is no such ambiguity for the local Chern marker, so that M C (r) is well defined at a specific r. Certainly the local conductivity tensor σ ij (r) = ∂j i (r)/∂E j , the first-order response of the local current density j(r) to a homogeneous macroscopic electric field E, is physically well defined. As shown by Rauch et al. [11], the local anomalous Hall conductivity, which is just the antisymmetric part of σ(r), can be written as a sum of the Chern marker of Eq. (15) and an additional non-geometric or "cross-gap" term. This analysis puts the definition of the local Chern marker M C (r) on a relatively firm footing.
Returning to the local markers for M LC and M IC , we impose four requirements, as follows.

Independence of origin
We first require that the trace projections should be independent of origin. A potential marker M [P,Q,H] (r) may be regarded as depending on the system-specific operators P , Q, and H as indicated by the superscript notation. Introducing the unitary operator T for a translation by displacement t, and defining translated oper-atorsP = T P T † and similarly for Q and H, the desired independence of origin is equivalent to asking that M [P,Q,H] (r) = M [P ,Q,H] (r + t). We accomplish this by insisting that the markers be written in terms of the operators X = P xQ, Y = P yQ, Z = P zQ (19) and their Hermitian conjugates. A combination like P rQ transforms toPrQ =P (r − t)Q =P rQ, where the last equality follows from the fact that P Q = 0. Thus, expressions built out of the operators in Eq. (19) will automatically satisfy the desired independence of origin. With this notation, Eqs. (13), (17), and (18) become Other choices such as r|yP xQHQ|r would not satisfy this property. At this point, then, we have three candidates for the LC marker, and there is a corresponding set of three choices for the IC marker.

Covariance under rotations
We next insist that the markers should transform as vectors under global rotations of the system. That is, given a spatial rotation R and its corresponding unitary operator R, we ask that In particular, we insist that any candidate for the marker z component should be invariant under rotation by an arbitrary angle θ aboutẑ. For θ = π/2 we have X → Y and Y → −X, and using the general formula nor M R LC satisfies rotational invariance by itself, but at least for θ = π/2, the average (M L LC + M R LC )/2 does. A more careful check confirms that this combination is also invariant under arbitrary rotations by angle θ. We can therefore propose two valid local-circulation markers, where the z subscript has temporarily been restored, and the E superscript in the second equation involving the anticommutator indicates that H is in an "external" (left or right) position. More generally, we can define vector markers M LC whose x and y components have the same form as in Eqs. (23) and (24) but with a cyclic permutation of the Cartesian indices. Since any rotation can be written as a product of three consecutive rotations around Cartesian axes by Euler angles, and since our markers have the correct rotation properties under each of these separately, it follows that the overall rotational properties are guaranteed. That is, the rotation of the marker is consistent with the rotation of the system. The same reasoning applies to the itinerant circulation markers where we have reverted to dropping the z subscript. Combining with Eqs. (23) and (24), we can arrive at several suitable expressions for the full marker in Eq. (12). If we consistently use the center or external versions, we arrive at two choices We also consider a third marker involving a mixed choice, We could also introduce its partner M E,C , but there are actually only three linearly independent markers, since it is easy to show that

Transformation under H → −H
Our third requirement is based on the following observation regarding the behavior of microscopic currents. As will be shown later in this section, the microscopic current density at r is given by Under a transformation that acts to change the overall sign of H, the occupied and unoccupied state manifolds are switched; that is, P = Q and Q = P . Substituting H = −H and P into Eq. (30), and using the facts that Q = 1 − P and that the expectation value of a Hermitian operator is real, it is easy to see that j(r) is left unchanged. Therefore, the surface magnetizations and hinge currents must be left unchanged as well. When H → −H, we see that X → X † , Y → Y † , and Z → Z † . It is then readily apparent that M C,C → M C,C and M E,E → M E,E . However, surface magnetization from the M C,E marker is not left invariant, as M C,E → M E,C . We therefore eliminate M C,E from the list of valid markers.

Magnetic quadrupole of finite systems
To conclude our list of acceptable marker properties, we turn to a discussion of the possible role of the bulk magnetic quadrupole moment (MQM) in a theory of surface magnetization. Our motivation for considering the MQM stems from the recently developed theory of boundary electric polarization [15,16]; for a crystallite with vanishing bulk polarization, the theory identifies the bulk electric quadrupole moment as contributing to the electric polarization at the surface. By analogy, it would then be unsurprising if the MQM similarly contributed to the magnetization at a surface for systems with zero bulk magnetization. Earlier works have focused on deriving expressions for the MQM within periodic boundary conditions [25][26][27], but have not discussed whether or how it manifests itself on the boundary of a bulk.
As a theory of the bulk MQM density for an extended system is not well established, we focus here on the case of an ideal molecular crystal, i.e., a collection of identical independent units, which we refer to as "molecules", arranged without overlap on a crystal lattice.
Classically, the MQM tensor Q ij of a single molecule is defined in terms of its microscopic current distribution j(r) as This is a traceless tensor, since Its antisymmetric part is known as the toroidal moment, while its symmetric part is solely responsible for the quadrupole contribution to the far B field in the multipole expansion. When this molecular unit is periodically repeated to generate a large but finite crystallite, the quadrupole moment densities are defined as Q = Q/V and Q = Q/V , where V is the unit cell volume. One can then show that the MQM generates a magnetization at a surface with outward normaln given by [28] Repackaging the toroidal part of the MQM density as the vector m i = ijk Q jk /2, it is clear that this produces a surface magnetization M surf =n × m that is parallel to the surface. This part does not contribute to the surfacenormal component M ⊥ , which is just given by We are thus free to insert either the unsymmetrized Q or the symmetrized Q tensor into Eq. (35) when computing M ⊥ values, ignoring the toroidal component [29,30]. We will refer to the MQM tensors defined as in Eq. (31) as current-based quadrupoles. In the special case that the current distribution takes the form of a collection of point magnetic dipoles m µ located at sites r µ , one can introduce an alternative dipole-based definition If one considers the limit that the point dipoles are constructed from small current loops of vanishing size, it turns out that the traceless part of Q dip is identical with Q (see Appendix C). That is, For a molecular crystal, the dipole formula provides a natural connection to the use of the local magnetization marker. That is, we can calculate the MQM density of the crystal as or where M i (r) and M µ,i are the continuum and discrete markers, respectively. Therefore, in the rest of this paper, we will focus on the dipole-based expression in Eq. (39) in the context of our tight-binding calculation of the local markers. We will then remove the trace part in order to compare with the current-based tensor computed from Eq. (31) in the same TB framework. Tests of this kind will be presented in Sec. V D. Unlike the other criteria described above in Secs. IV A 1 to IV A 3, we do not know how to determine a priori whether a certain marker will satisfy the present criterion of reproducing the MQM correctly for molecular crystals. However, we will see in Sec. V D that some choices of marker sometimes fail this test. Thus, even if numerical in nature, these tests provide important constraints on the appropriate choice of marker for the computation of surface magnetization.
We note in passing that the use of the raw dipole-based MQM, without removing the trace part, will generate different results for the surface magnetizations M ⊥ . However, it will shift M ⊥ equally on all facets, and therefore will still yield a correct prediction of the hinge currents, as these only depend on differences of M ⊥ on neighboring facets (see Eq. (2)). In this sense, the use of the raw dipole-based MQM in place of the traceless one is reminiscent of the "shift freedom" in the definition of M ⊥ discussed in the Introduction, here for the case of an idealized molecular crystal.

Discussion
In later sections we find it convenient to work with the two markers given in Eqs. (27)(28) above. The marker M C,C of Eq. (27) is the original one of Bianco and Resta [9], while M E,E is a new candidate that has not, to our knowledge, been considered previously.
To end the discussion of the markers, we also note that there is no restriction on using a linear combination of M C,C and M E,E to write down a new marker, as long the linear combination reproduces the bulk magnetization M of Eq. (8). To maintain this requirement, the (real) coefficients a CC and a EE of the linear combination must sum to unity.
We select the particular linear combination (a CC , a EE ) = (1/3, 2/3) to form the "linear" marker This marker is also expressible as Having formulated our list of markers, we are faced with several questions regarding their behavior. First, it remains to be seen which, if any, of the markers will correctly predict the hinge current. For each surface forming a hinge, identical markers will be used to compute their respective magnetizations, which in turn will be substituted into Eq. (2) to yield the current. If the system only exhibits axion-even symmetries, and the markers disagree on the values of the surface magnetizations, but their differences yield the hinge current correctly, then this would be a realization of the shift freedom discussed earlier. Second, if the system does feature an axion-odd symmetry, will all three markers yield the same value of the surface magnetization? If not, then which of the markers correctly predict the hinge current?
Before attempting to answer these questions, we discuss some technicalities of the computation of the surface magnetization from the local markers.

B. Surface magnetization and macroscopic averaging
Recall that we assume that the system has enough symmetry to force the macroscopic magnetization to vanish in the bulk. However, this does not necessarily imply that M(r) vanishes everywhere, only that its unit cell average vanishes. In case M(r) does vanish identically deep in the bulk, then for any given trace projection, it is straightforward to integrate M(r) over one surface unit cell, with the surface-normal integral carried deep enough into bulk and vacuum, to obtain a corresponding macroscopic surface magnetization. If this is not the case, however, a coarse-graining procedure is required, as described next.
We construct an insulating slab of thickness L with the outward normals at the top and bottom surfaces beinĝ z and −ẑ respectively, and with cell-periodic boundary conditions in x and y. The slab features the maximal symmetry allowed by the bulk. We assume L >> c, where c is the lattice constant in the surface-normal direction. We define a "layer-resolved" local marker by averaging the local marker over an in-plane unit cell of area A at fixed z: Note that the total magnetization of the slab (magnetic moment per unit area) is the integral of M(z) over the thickness of the slab. The macroscopic surface magnetization is determined from the application of a smoothening procedure for the layer-resolved marker, since simple sums of the marker may not be convergent. We employ the sliding window average approach to compute the surface magnetization; for details on this method, we refer the reader to Sec. IV.A of Ref. [16]. In numerical work, the averaging amounts to an integration of the local marker weighted by a "ramp function" that extends sufficiently deep in the bulk. The ramp-down function is defined as Letting the bottom surface of the slab be located at z = 0, we then get where the range of z integration runs from deep in the vacuum to deep in the interior of the crystal. The result is independent of z 0 so long as it is sufficiently deep in the bulk region of the slab. The formalism for the calculation of a surface magnetization has been presented here in the context of a continuum-space treatment. However, our numerical tests will be performed on discrete TB models; in this setting, Eqs. (42) and (44) are expressed as discrete sums over TB sites, where the local markers are now defined. For further details, we refer the reader to Appendix A.

C. Macroscopic hinge current and macroscopic averaging
A direct calculation of the macroscopic hinge current I is necessary in order to test whether it is correctly pre- dicted by each of the various trace projections. The calculations will be accomplished via a macroscopic averaging procedure analogous to that of Sec. IV B, but now applied to the integration of the microscopic current density over the appropriate hinge region. In this section and in the rest of the paper, we label vector operators with a hat in order to distinguish them from ordinary vectors.
The microscopic current density j(r) is obtained from the expectation value of the microscopic current density operatorĵ(r). Given a single-particle Hamiltonian H, this is given byĵ Our calculations of the macroscopic current are performed on a series of insulating infinite rod geometries that retain periodic boundary conditions along one of the three Cartesian directions, as illustrated in Fig. 7. For concreteness, in this subsection we focus on computing the hinge current I FB for a rod running parallel to the x-axis, with unit-cell periodicity a. Let the rod have width W along y and height L along z, with W and L much larger than the cell dimensions b and c along these respective dimensions.
In analogy with the surface magnetization, the macroscopic current is obtained from a sliding window average performed over the microscopic currents at the corresponding hinge, but now in two dimensions. The appropriate weight function w(y, z) is a product of two ramp functions in y and z, namely w(y, z) = f b (y−y 0 )f c (z−z 0 ) where y 0 and z 0 are located sufficiently deep within the interior of the rod. I FB is then given by where the range of integration is over all y and z. The currents I BL and I LF are computed analogously with the appropriate cell parameters and Cartesian variables used in Eq. (47). As in the previous subsection, the formalism here is developed in the context of a continuum. In a TB formulation, Eqs. (46) and (47) will turn into discrete sums over TB sites. We refer the reader to Appendix B for further details.

V. NUMERICAL RESULTS
To test the results of the symmetry analysis of Sec. III and the local marker calculation of surface magnetization, we study three tight-binding (TB) models of spinless electrons. The first model is comprised of an infinite stack of the Haldane-model layers, while the other two models are composed of stacks of two-dimensional (2D) square-plaquette layers. Each model features symmetryenforced zero bulk magnetization and is considered at half-filling. The TB basis orbitals are assumed to have no orbital moment of their own and to diagonalize the position operators.
We focus on the variation of a single parameter in the Hamiltonian away from the special value that quantizes the CSA coupling, and investigate how the surface magnetization found from the various trace projections behaves as a function of this variation. For each model, the Fermi energy is set to zero, allowing us to neglect the Chern marker contribution to the local magnetization marker.
We use the slab and rod geometries outlined in Secs. IV B and IV C to compute the surface magnetizations and hinge currents, respectively. Fig. 7 illustrates the particular magnetizations and hinge currents that we compute. The surface magnetizations are found from the coarse-graining procedure described in Sec. IV B while the macroscopic hinge currents are found from the coarsegraining procedure of Sec IV C. Slab and rod geometries for the computation of surface magnetizations and hinge currents are generated by truncating the bulk, i.e., the hoppings to vacant sites are removed while other hoppings and site energies are unchanged. The electronic Hamiltonians for bulk, slab, and rod geometries are constructed and solved using the PythTB code package [31].
Details of our numerical results will be given in Appendix D. For each model, however, we provide here two tables to demonstrate the behavior of the surface magnetization and hinge currents for the axion-odd regime, as well as a representative axion-even setup. Each appropriate directly computed hinge currents displayed in the bottom-most rows of the tables. At the end of this section, we numerically address the question of which local marker is able to reproduce the current-based MQM for a finite system with no magnetic moment. We perform a test on a finite cubic model, and compare the current-based MQM tensor to the dipolebased MQM tensors derived from the M C,C , M E,E , and M lin local markers.

A. Alternating Haldane model
The first model we study consist of half-filled Haldanemodel layers [32] placed directly on top of each other, which are then subsequently coupled via interlayer hoppings along the (001) direction. The interlayer hoppings couple sites on identical sublattices, and alternate in value from one interlayer region to the next. The intralayer parameters are chosen to be such that if the layers were decoupled, their Chern numbers would vanish and their orbital magnetizations would form an up-down   The symmetries of this model ensure that the bulk magnetization is zero. In particular, there is a three-fold rotation axis C 3z about the lattice sites (or honeycomb centers) as well as a time-reversed mirror plane m z passing between the Haldane layers that eliminate the magnetization. Additionally, a time-reversed mirror plane m x bisects the honeycombs along their main diagonals, and a two-fold rotation axis C 2y lies between the layers and in the plane of the m x mirror.
The quantization of the CSA coupling is controlled by varying t 3 while keeping all other parameters fixed. When t 3 = t 3 , the system gains a time-reversed halftranslation {E |c/2} symmetry as well as an m z mirror plane residing in the layers; both symmetries are axionodd.
We The macroscopic hinge current I FB is found from an x-extensive rod composed of twenty cells along a 3 and twenty cells along a 2 . The y-extensive rod employed for I BL is composed of twenty cells along both A 1 and A 3 , while the z-extensive rod for I LF is composed of twenty cells each along A 1 and A 2 . For each rod, a 1D k-space sampling of thirty points is employed.
For t 3 = 0.1 and t 3 = 0.15, respectively, Tables III and IV report the values of the surface magnetizations computed from the different markers. From the tables, we note that for a given surface all three markers yield identical magnetizations, independent of whether or not  the system features an axion-odd symmetry. Furthermore, the magnetizations always yield the correct values of the hinge currents. In the axion-odd case of t 3 = 0.1, we observe that M F ⊥ and M L ⊥ may be understood to be zero due to the slab constructions for their calculations featuring the m z mirror plane. By the same token, I LF vanishes as well, as the z extensive rod used for its calculation also features the same mirror symmetry.

B. Two-layer square-plaquette model
The model studied here draws inspiration from a 2D square-plaquette model developed in Ref. [5]. As presented there, the model consists of a nearest-neighbor TB Hamiltonian whose primitive cell consists of four sites labeled A-D as shown in Fig. 9(a), with lattice vectors a 1 =x and a 2 =ŷ. Each of the four resulting plaquettes in the unit cell is threaded by a magnetic flux that corresponds to endowing some of the hopping amplitudes with an identical complex phase e iφ/2 . The moduli of all the nonzero hoppings are set to an identical value t. Fig. 9(b) provides a detailed illustration of the model that we study. Our model features equidistantly spaced layers of the plaquette model stacked directly on top of one another along theẑ direction. Subsequent layers are  This model possesses a time-reversed inversion symmetry I about a point located midway between C sites on neighboring layers, as well as a two-fold rotation C 2xy about an axis located between layers and passing above A and C sites. Additionally, there are m xy and m xy mirror planes that pass through sites B and D, and A and C, respectively. There is also a time reversed glide mirror plane {m z |τ /2} passing in between layers with τ =x+ŷ. These symmetries ensure that the bulk magnetization remains zero.
The CSA coupling is tuned by keeping all parameters fixed except t 3 . If t 3 and t 3 are equal, the model gains an axion-odd {C 2z |c/2} screw axis passing through the C sites. The parameters for this model are chosen as ∆ = −0.5, t = 0.06, φ = π/4, t 3 = 0.1, and 0.05 ≤ t 3 ≤ 0.15.
We will consider two distinct surface terminations in the directions transverse to the layer stacking, and for each we will investigate the behavior of the resulting surface magnetizations and hinge currents. The first surface termination results from creating slabs along a 1 and a 2 , and we refer to these as the straight-edge terminations. The second type of termination results from cutting along the a 1 ±a 2 directions, which we refer to as the zigzag termination. For our calculations with this termination, we construct supercells with lattice vectors A 1 = a 1 − a 2 , A 2 = a 1 + a 2 , and A 3 = a 3 , and cut slabs along the A 1 and A 2 directions. In this case M B ⊥ , M F ⊥ , and M L ⊥ are found from slabs with eleven unit cells along a 3 , seven unit cells along a 2 , and seven unit cells along a 1 , respectively. For each magnetization, a 7 × 7 k-space mesh in reduced coordinates is used.
The macroscopic hinge current I FB is found using an xextensive rod composed of eleven unit cells along a 3 and seven unit cells along a 2 . The y-extensive rod employed for I BL is composed of eleven unit cells along a 3 and seven unit cells along a 1 , while the z-extensive rod for I LF is composed of seven unit cells along both a 1 and a 2 . For each rod, a 1D k-space sampling of thirty points is employed.
For In the axion-odd regime when t 3 = 0.1, all markers yield the same value of magnetization at a single surface. However, in the axion-even regime when t 3 = 0.15, the surface magnetizations derived from different markers do not agree, and furthermore only the M lin marker correctly predicts the hinge currents. We will soon see that among the three markers, it appears that only M lin consistently predicts the hinge currents, regardless of whether the system is in the axion-odd or axion-even regime.  In order to describe the slabs and rods that occur when zigzag terminations are present, we turn to using the supercell vectors A 1 , A 2 , and A 3 described at the beginning of the subsection. We align the x and y axes along A 1 and A 2 respectively. Since the bottom surface of a slab cut along A 3 is identical up to a rotation in the x-y plane to the bottom surface of the slabs used to compute M B ⊥ in the straight-edge case, the values of M B ⊥ will be left unchanged. With the rotation of the coordinate axes, we also observe that the m xȳ and m xy mirrors become m x and m y mirrors, respectively, the C 2xȳ axis becomes a C 2x axis, and {m z |τ /2} becomes {m z |A 2 /2}.
In this case M F ⊥ and M L ⊥ are found from slabs with four cells along A 2 and A 1 , respectively. For each magnetization, a 7 × 7 k-space mesh in reduced coordinates is used. The macroscopic hinge current I FB is found from an xextensive rod composed of nine cells along A 3 and nine cells along A 2 . The y-extensive rod employed for I BL is composed of nine cells along A 3 and nine cells along A 1 , while the z-extensive rod for I LF is composed of nine and nine cells along A 1 and A 2 , respectively. For each rod, a 1D k-space sampling of five points is employed.
Tables VII and VIII demonstrate the values of the surface magnetizations computed from the different markers for t 3 = 0.1 and t 3 = 0.15, respectively. In the axion-odd case, we see that for a fixed surface all the markers yield identical values of the magnetization. Furthermore, the values of M F ⊥ and M L ⊥ vanish due to the m x , m y , and {C 2z |c/2} symmetries. m x ensures identical M ⊥ values on surfaces with outward unit normalsx and −x, while m y does so for surfaces with outward unit normalsŷ and −ŷ. {C 2z |c/2} then ensures that these surface magnetizations are zero.
In the axion-even regime, the markers at a single surface facet yield differing values of the magnetizations, and none except the 'lin' marker predict the hinge currents correctly.
These observations are highlighted in Fig. 10, which displays the hinge currents as a function of t 3 and compares the currents predicted using the different markers to the directly computed hinge currents. The model studied in this section employs the same underlying square-plaquette layers that were used in Sec. V B and is depicted in Fig. 9(c). Just like the model of Sec. V B, the present model features equidistantly spaced layers of the plaquette model stacked directly on top of one another along theẑ direction, but now the adjacent layers are related by a time-reversed rotation C 4z about a z axis passing through the C site, resulting in four layers per unit cell. The lattice vector along the stacking direction is a 3 =ẑ. Otherwise, the models share the same features.
This model features a two-fold rotation axis C 2x between the second and third layers in the unit cell and passing above the C and D sites. There is also a {C 2z |c/2} screw axis that passes through the C sites. Together these symmetries eliminate the bulk magnetization. If the interlayer hoppings t 3 and t 3 are equal, the model gains an additional axion-odd {C 4z |c/4} screw axis running through the C sites. All parameters except t 3 are kept fixed in order to access the quantized CSA coupling  For this model, we also study the behavior of hinge currents and surface magnetizations for the different transverse surface terminations.

Four-layer plaquette model: straight edges
The magnetizations M B ⊥ , M F ⊥ , and M L ⊥ are found from slabs with six unit cells along a 3 , four unit cells along a 2 , and four unit cells along a 1 , respectively. For each magnetization, a 7 × 7 k-space mesh in reduced coordinates is used.
The macroscopic hinge current I FB is found using an xextensive rod composed of six unit cells along a 3 and four unit cells along a 2 . The y-extensive rod employed for I BL is composed of six unit cells along a 3 and four unit cells along a 1 , while the z-extensive rod for I LF is composed of four unit cells along both a 1 and a 2 . For each rod, a 1D k-space sampling of thirty points is employed.
For t 3 = 0.1 and t 3 = 0.15, respectively, Tables IX and X report the values of the surface magnetizations The same is true in the axion-even regime of Table X. Fig. 11 demonstrates these facts for the prediction of the I LF hinge current. For all the values of t 3 that we used, covering both axion-odd and axion-even regimes, only the 'lin' marker was consistent in its prediction of the current.

Four-layer plaquette model: zigzag edges
As in the case of two-layer plaquette model with zigzag terminations in Sec. V B 2, we turn to using the supercell  vectors A 1 = a 1 − a 2 , A 2 = a 1 + a 2 , and A 3 = a 3 for our calculations involving the zigzag edges of the four-layer plaquette model. We align the x and y axes along A 1 and A 2 respectively. Since the bottom surface of a slab cut along A 3 is identical up to a rotation in the x-y plane to the bottom surface of the slabs used to compute M B ⊥ in the straight-edge case, the values of M B ⊥ will be left unchanged. With the rotation of the coordinate axes, C 2x becomes C 2xy . M F ⊥ and M L ⊥ in this case are found from slabs with three cells along A 2 and A 1 , respectively. For each magnetization, a 5 × 5 k-space mesh is used in reduced coordinates.
The macroscopic hinge current I FB is found from an x-extensive rod composed of four unit cells along A 3 and six unit cells along A 2 . The y-extensive rod employed for I BL is composed of four unit cells along A 3 and six unit cells along A 1 , while the z-extensive rod for I LF is composed of six cells along both A 1 and A 2 . For each rod, a 1D k-space sampling of 5 points is employed.
Tables XI and XII report the values of the surface magnetizations computed from the different markers for t 3 = 0.1 and t 3 = 0.15, respectively. This time, in the axion-odd regime when t 3 = 0.1, we see that for any given surface facet, all markers yield the same value of magnetization, and when compared across facets, all accurately predict the hinge currents. The values of M F ⊥ and M L ⊥ vanish, which may be understood to be a consequence of the simultaneous C 2xy and {C 4z |c/4} symmetries. The former symmetry implies that M F ⊥ = M L ⊥ , while the latter results in both vanishing.
In the axion-even regime when t 3 = 0.15, we find a result that is more reminiscent of the straight-edge case. That is, even for a single surface facet, the markers yield different values of the magnetization, and the hinge currents computed from the differences of M ⊥ values disagree with each other when computed using different markers. They also disagree with the direct calculation of the hinge current except for one case, namely that of the M lin marker. M F ⊥ and M L ⊥ are found to be identical, as the C 2xy symmetry is still present in this setup. However, only the M lin marker yields magnetizations that correctly predict the hinge currents.

Summary
For the four-layer plaquette model, unlike the models studied previously, we find that different markers generally disagree on the values of M ⊥ and on their predictions for the hinge currents. This occurs in both the axion-odd and axion-even cases for straight-edge surface terminations, but only in the axion-even case for zigzag surface terminations. Only the M lin marker consistently predicts the correct hinge currents.

D. Local markers and magnetic quadrupole moment of a finite system
We now return to investigate the questions posed in Sec. IV A 4 by conducting a numerical study of a spinless TB model of a cubic "molecule". The system is designed to have a symmetry that enforces a vanishing magnetic dipole moment, but otherwise is as arbitrary as possible. In particular, it has no axion-odd symmetries. The structure, TB site labeling, and onsite energies of the model are depicted in Fig. 12. Additionally, the model features hoppings t ij to nearest, second-nearest, and third-nearest neighbors, where the subscripts indicate a hopping from TB site i to TB site j. The hoppings are generally complex with t ij = t * ji to ensure that the Hamiltonian is Hermitian.
The Hamiltonian is chosen to obey an S 4z = m z C 4z improper rotation symmetry, where the m z mirror plane bisects the molecule, while the C 4z axis runs through the center of the cube. The resulting magnetic point group of the system is then {E, S 4z , (S 4z ) −1 , C 2z }, where E is the identity operation and C 2z is the square of S 4z . These symmetries ensure a vanishing magnetic dipole moment, resulting in an origin-independent, physically meaningful MQM for the system. The magnetic point group also ensures that the MQM features Q xx = Q yy and Q xy = −Q yx . Additionally, Q zz = −2Q xx , as the tensor is traceless.
In addition to the onsite energies, the Hamiltonian is parametrized by eight hopping parameters. The nearest neighbor hopping terms t 1 , t 2 , and t 3 parametrize, respectively, the hoppings

VI. DISCUSSION
Based on the symmetry analysis of Sec. III, we expect that surface magnetization should be well defined only when a system exhibits an axion-odd symmetry. More generally, we expect that differences of M ⊥ values on neighboring facets should correctly predict the bound current flowing on the hinge where the facets meet, even when axion-odd symmetries are absent.
Focusing first on the axion-odd case, we find that both expectations are sometimes violated by some choices of marker. We found no such violation for the alternating Haldane and two-layer square-plaquette models. However, the four-layer square-plaquette model does show signs of trouble. For the straight-edge surface terminations, the magnetizations found from the M C,C , M E,E , and M lin markers do not agree with each other. Furthermore, among these, we found that only the M lin marker correctly predicts the hinge current. Thus, while the physical magnetization is uniquely defined in the axionodd case, only the M lin marker consistently gives the correct value for it in our tests.
When tuning from the axion-odd to the axion-even regime, the same conclusion was reinforced. While all markers continued to agree with each other for the stacked Haldane model, they disagreed in all other cases, i.e., for the two-layer square-plaquette model and all surface terminations of the four-layer square-plaquette model. Again, the only marker that always reproduced the hinge currents across all models was the M lin marker.
The importance of the M lin marker was also highlighted by our numerical finding that the dipole-based traceless MQM derived from it appears to correctly reproduce the current-based MQM in finite systems. Thus, we provisionally identify M lin as the best candidate for a marker-based formulation of a theory of surface orbital magnetization.
Given the M lin marker's evident connection to the MQM, we have also explored the symmetries of our models to understand whether they permit a nonzero bulk MQM. While the question of defining an orbital MQM for a general bulk material has received some attention [25,26], the situation appears to remain unsettled. Assuming such a definition is possible, however, it must respect the symmetries of the magnetic point group.
The four-layer square-plaquette model is most informative in this respect. Focusing first on the axion-odd case, the tensor is diagonal with Q xx = −Q yy nonzero and Q zz = 0 when the coordinate axes are chosen as in Fig. 9(c) for the straight-edge surfaces. In the zigzagedge case, we rotated the axes aboutẑ by −π/4, transforming the tensor such that the only nonzero elements were Q xy = Q yx . We found that the markers disagreed with each other in precisely those cases in which the corresponding diagonal MQM tensor element was nonzero, namely for the side surfaces in the straight-edge case. According to Eq. (34), these were also the cases in which the bulk MQM has the right symmetry to contribute to the surface magnetization. When the axion-odd symmetry was broken, we found that all three diagonal elements were present in the MQM for both surface terminations, and that the markers disagreed with each other for all surfaces For the alternating Haldane and two-layer squareplaquette models, the symmetry of the MQM is such that it vanishes in the axion-odd regime but acquires some diagonal elements in the axion-even regime. For the latter model, the markers also failed to agree in the axion-even regime.
To summarize this discussion, when inspecting all three models, we found that the only cases in which markers disagree, with only the M lin correctly predicting the hinge currents, were those in which diagonal elements of the MQM tensor were present. Thus, however the bulk MQM might be defined, our results suggest that only the M lin marker accurately takes into account its contribution to surface magnetization.
Notably, our findings provide no evidence for the possible marker "shift freedom" discussed in Sec. I. That is, we did not find any two markers that always yielded the correct hinge currents, but resulted in surface magnetizations differing from each other by a constant shift in axion-even cases.
We find ourselves, then, in a somewhat ambiguous position. Our empirical evidence, based on the model studies, indicates that the original Bianco-Resta marker M CC is not suitable for coarse-graining to obtain the surface magnetization. Instead, a modification of it, namely M lin , has passed all tests and appears to be a suitable marker. Still, despite a concerted effort to find a formal connection, we do not yet have a fundamental understanding as to why M lin succeeds where the others do not. Our work thus raises this issue as being of paramount importance to be addressed by future theoretical investigations.
Of course, it would be desirable to make contact with experiment. We have argued that the hinge currents are observable in principle via the magnetic fields they generate, which could be observed by scanning magnetic probes. The surface magnetization M ⊥ itself in the interior of a facet does not produce any observable electric or magnetic potential, and so cannot be observed directly. It is possible that optical probes, such as terahertz Kerr reflectivity, might provide information. However, one should keep in mind that when symmetry allows for the presence of a surface orbital magnetization, it also allows for a surface spin magnetization, which is likely to be much larger. Thus, filtering out just the orbital component may be challenging.
Other generalizations of our present work remain to be developed. The symmetry analysis of Sec. III was con-ducted in full generality and is applicable to insulating as well as metallic systems at any temperature. However, it is not immediately clear whether it is possible to compute surface magnetization using the local marker for T = 0. A crucial assumption going into the construction of the local marker is the idempotency of P [9]; in the case of finite temperature, P no longer displays this behavior. Additionally, for metals at T = 0, even though P remains idempotent, it is no longer exponentially localized, with P (r, r ) featuring power law decay with |r − r |. It is not immediately clear how this affects the local marker and its subsequent coarse-grained average. In these situations, it appears likely that a different method of computing surface magnetization will be needed.
Additional directions to explore include the presence of facet magnetizations for systems with bulk magnetization, as well as for axion-odd systems with θ CS = π. It would also be interesting to see whether a Wannier-based formulation for surface magnetization is possible, at least in insulators. This is particularly attractive as the theories of bulk electric polarization, bulk orbital magnetization, and edge electric polarization can all be developed using a Wannier-based approach [1,2,4,5,16,17]. Finally, we mention the prospects for a first-principles density-functional theory (DFT) implementation of the calculation of surface magnetization for topologically trivial bulk materials. The output of any DFT calculation may be adapted to a TB framework via Wannier interpolation, which may be accomplished by code packages such as Wannier90 [33]. It is important to keep in mind, however, that in our calculations of the local markers we employed the diagonal approximation for the position operator. In the case of a Wannier interpolation, the position operator is not necessarily diagonal in the basis of Wannier states, and the off-diagonal terms must be kept in mind when performing calculations.

VII. SUMMARY
In this work, we have explored the possibility of defining surface orbital magnetization for insulating systems that are topologically trivial in the bulk and feature no bulk orbital magnetization. We have demonstrated that the physical observable corresponding to surface orbital magnetization is the current running along a hinge shared by two surface facets. As a result, only differences between magnetizations at adjacent surface facets are well defined in general. By means of a symmetry analysis, we have shown that individual values of facet magnetizations should be well defined when the bulk symmetry group contains an axion-odd symmetry. We further expanded on the local marker formulation of orbital magnetization developed by Bianco and Resta, and introduced a set of markers satisfying a list of physically meaningful transformation requirements. The coarse-grained averages of these markers were then used to compute surface magnetizations.
We tested the conclusions of our symmetry analysis and our formalism for computing surface magnetization on a series of spinless TB models. We found that the markers do not always agree on the value of the surface magnetization, even in the axion-odd case where we know from symmetry arguments that M ⊥ is well defined. Our results indicate that only a single marker on our list, the M lin marker, consistently predicts the correct hinge currents. According to the symmetry considerations of Sec. III, in the axion-odd case we can then conclude that it has correctly computed the unique M ⊥ values. In the axion-even case, instead, we can no longer say that it predicts uniquely correct M ⊥ values, only correct differences of M ⊥ values.
We additionally tested the markers to understand whether any of them correctly generated the currentbased MQM of Eq. (31), which we argued may contribute to surface magnetization. We found that only the M lin marker correctly reproduced the current-based MQM when applied to a finite system with a nontrivial MQM and vanishing magnetic dipole moment. We also observed that when the different markers failed to agree on the surface magnetizations for axion-odd systems, the bulk symmetry group permitted an MQM that could contribute to the surface magnetizations.
Overall, our study indicates that a coarse-graining of the M lin marker provides a suitable framework for computing surface magnetizations. However, a formal derivation that would explain the unique success of the M lin marker remains elusive, and questions persist about the definition and role of the bulk quadrupole tensor. Thus, while our work establishes a foundation for a theory of surface orbital magnetization, it also highlights the need for further work to resolve some important remaining questions.
This work was supported by NSF Grants DMR-1954856 and DGE-1842213. The calculations for all other LC and IC markers are similar, and we will simply state the results at the end of this Appendix. We remind the reader that in all of our models, we have set the Fermi energy inside the band gap at zero, so that we ignore the Chern marker contribution to the local marker. For details on how to compute the Chern marker in the TB representation, we refer the reader to Ref. [20].
We denote the valence (v) and conduction (c) band eigenstates of the slab Hamiltonian as ψ vk (r) and ψ ck (r), respectively, where k = (k x , k y ). The valence and conduction band projectors are then written as P = (1/N k ) vk |ψ vk ψ vk | and Q = (1/N k ) ck |ψ ck ψ ck |, where N k is the number of k points in the 2D Brillouin zone mesh. With this in mind, we observe that where E nk is the eigenenergy corresponding to |ψ nk , X vck = ψ vk |x|ψ ck , and Y † cvk = ψ ck |y|ψ vk . In our calculations of the local markers, we do not directly compute the matrix elements of the position operators x and y in the energy eigenbasis; rather, we use the velocity operatorv, defined aŝ to rewrite and similarly for Y † cv k . The three LC markers are then The three IC markers are The combinations of these markers then yield the M C,C , M E,E , and M lin markers.
elements of the position operators in the TB basis, where r i denotes the position of TB orbital i. With this approximation in mind, the matrix elements of the velocity operatorv = (1/i )[r, H] in the TB basis are expressed as Given a set of occupied states |ψ n , the total current in the system is given by J = describes the total current flowing on the bond ij . Since the current is assumed to flow uniformly along the straight-line bond connecting the sites, it makes no difference if we treat it as distributed uniformly along the length, concentrated at the center, or partitioned between the two end points of the bond. We make the last choice in order to write the total current as a sum over sites, J = i j i , where We treat this as the discrete analog of Eq. (46), and we use it in the calculations of the macroscopic hinge currents in our TB models by performing the coarse-graining procedure of Sec. IV C on the j i .
Finally, we note that in the case of the extended rod geometries used to find the macroscopic hinge currents in the main text, the valence projector P is written in terms of the rod valence eigenstates ψ vk as P = (1/N k ) vk |ψ vk ψ vk |, where the sum is over the N k points of the mesh spanning the 1D Brillouin zone. This expression for P is then substituted into Eq. (B6).