Collinear Orbital Antiferromagnetic Order and Magnetoelectricity in Quasi-2D Itinerant-Electron Paramagnets, Ferromagnets and Antiferromagnets

We develop a comprehensive theory for magnetoelectricity in magnetically ordered quasi-2D systems whereby in thermal equilibrium an electric field can induce a magnetization $m$ and a magnetic field can induce a polarization. This effect requires that both space-inversion and time-reversal symmetry are broken. Antiferromagnetic (AFM) order plays a central role in this theory. We define a N\'eel operator $\tau$ such that a nonzero expectation value $\langle \tau \rangle$ signals AFM order, in the same way $m$ signals ferromagnetic (FM) order. While $m$ is even under space inversion and odd under time reversal, $\tau$ describes a toroidal moment that is odd under both symmetries. Thus $m$ and $\langle \tau \rangle$ quantify complementary aspects of magnetic order in solids. In quasi-2D systems FM order can be attributed to dipolar equilibrium currents that give rise to $m$. In the same way, AFM order arises from quadrupolar currents that generate the moment $\langle \tau \rangle$. The electric-field-induced magnetization can then be attributed to the electric manipulation of the quadrupolar currents. We develop a $k \cdot p$ envelope-function theory for AFM diamond structures that allows us to derive explicit expressions for the operator $\tau$. Considering FM zincblende and AFM diamond, we derive quantitative expressions for the magnetoelectric responses due to electric and magnetic fields that reveal explicitly the inherent duality of these responses required by thermodynamics. Magnetoelectricity is found to be small in realistic calculations for quasi-2D electron systems. The magnetoelectric response of quasi-2D hole systems turns out to be sizable, however, with moderate electric fields being able to induce a magnetic moment of one Bohr magneton per charge carrier. Our theory provides a broad framework for the manipulation of magnetic order by means of external fields.


I. INTRODUCTION
The technological viability of alternative spin-based electronics prototypes [1][2][3] hinges on the ability to efficiently manipulate magnetizations using electric currents or voltages. Various basic device architectures are currently being explored that could offer the crucially needed electric magnetization control. One promising approach utilizes antiferromagnetic materials [4,5], while another employs spin torques generated via the Edelstein effect [6][7][8][9]. A third interesting avenue has been opened by harnessing the magnetoelectric effect [10][11][12][13][14] in multiferroic materials [15][16][17][18][19] for switching the magnetization of an adjacent ferromagnetic contact [20,21]. Results obtained in our work point to an appealing alternative possibility, whereby intrinsic magnetoelectric couplings in ferromagnetic and antiferromagnetic quasi-twodimensional (quasi-2D) itinerant electron systems provide a nondissipative mechanism for electric control of magnetizations. * Dedicated to Ulrich Rössler on the occasion of his 80th birthday.
We present a comprehensive theoretical study of magnetoelectricity in these paradigmatic nanoelectronic structures that have the potential to become blueprints for future spintronic devices.
The magnetoelectric effect [10][11][12][13][14] refers to the property of certain materials to develop an equilibrium magnetic response to an electric stimulus, and vice versa. A systematic understanding of magnetoelectric phenomena can be based on an expansion of the free-energy density F as a function of the externally applied electric field E and magnetic field B [11,14], The first two lines in Eq. (1) pertain to ordinary electromagnetic phenomena [22], whereas terms in the third line are associated with magnetoelectricity. In particular, the magnetoelectric tensor α ij characterizes the generation of an electric polarization by a magnetic field and of a magnetization by an electric field, as is clear from the explicit expressions for the polarization P = −∂F/∂ E, . . , (2a) and the magnetization M = −∂F/∂B, Here and in the following, we have denoted by ∂/∂a the gradient vector (∂ ax , ∂ ay , ∂ az ) of derivatives w.r.t. the Cartesian components of a vector a ≡ (a x , a y , a z ). In both Eqs. (2a) and (2b), the first line embodies conventional electromagnetism in the solid state [11], whereas terms in the second line of these equations are ramifications of the magnetoelectric effect [10,11]. The appearance of the same set of coefficients α ij , β ijk , and γ ijk in these equations indicates a deep connection between the microscopic mechanisms causing a magnetically induced polarization and the microscopic mechanisms causing an electrically induced magnetization. As shown in the present work, quasi-2D systems facilitate the detailed discussion and thorough elucidation of the underlying mechanisms for such dual magnetoelectric responses. They also present a promising platform for exploiting magnetoelectricity in device applications.
As the product of E and B is odd under space inversion and time reversal, a nonzero tensor α ij is permitted only for systems with space-inversion symmetry and time-reversal symmetry both broken [11]. Terms proportional to the tensors β ijk and γ ijk embody higher-order magnetoelectric effects [14,23,24]. Systems in which only space-inversion (time-reversal) symmetry is broken can have nonzero tensors β ijk (γ ijk ), while α ij = 0. As an example for the latter in the context of the present work, we show that paramagnetic quantum wells in zincblendestructure materials exhibit the higher-order magnetoelectric effect associated with the tensor β ijk .
The magnetoelectric effect has been studied experimentally for a range of materials including ferromagnetic, antiferromagnetic and multiferroic systems [10,14,15,25,26]. Existing theoretical studies of the magnetoelectric effect have either focused on elucidating general properties of the tensors α ij , β ijk and γ ijk based on symmetry [27][28][29] or developed first-principles methods for their numerical calculation [30][31][32][33][34][35]. Typically, these works have also limited their scope to investigating only one of the two dual magnetoelectric responses. As a result, the microscopic basis for the intrinsic symmetry of electric and magnetic responses has been rarely discussed [36]. In contrast, the conceptually transparent and practically important quantum-well system considered in the present work provides a versatile, unified theoretical framework for describing magnetoelectricity in paramagnets, ferromagnets and antiferromagnets, covering both the electrically induced magnetization and the magnetically induced polarization and demonstrating explicitly how these two effects are intrinsically related.
Our realistic theoretical study focuses on the technologically important class of materials realizing variants of the diamond structure; see Fig. 1. As discussed earlier, magnetoelectricity only occurs in situations where both space-inversion and time-reversal symmetry are broken. Hence, the magnetoelectric effect is absent in paramagnetic materials having the inversion-symmetric [37] diamond structure [ Fig. 1(a)]. In contrast, the zincblende structure [ Fig. 1(b)] breaks inversion symmetry. In addition, time-reversal symmetry is broken in magnetized samples with ordered spin magnetic moments or with an orbital magnetization due to dissipationless equilibrium currents. Such a magnetization can be caused by a Zeeman coupling of the charge carriers to an applied magnetic field, or by a ferromagnetic exchange field [38,39] that is present in the material itself or induced by proximity to a ferromagnet. The origin of the magnetization is largely irrelevant for the microscopic mechanism of magnetoelectricity so that we denote all these scenarios jointly as ferromagnetically ordered. We demonstrate in this work the emergence of finite magnetoelectric couplings in ferromagnetically ordered quantum wells made from materials having a zincblende structure. We find that already in the absence of external fields, the interplay of broken space-inversion and time-reversal symmetry generates a collinear orbital antiferromagnetic order of the charge carriers that renders these systems to be actually ferrimagnetic. The magnetoelectric effect can then be viewed as arising from the manipulation of the equilibrium current distributions underlying the orbital antiferromagnetic order. Specifically, an electric field affects these currents in a way reminiscent of the Lorentz force such that the modified currents give rise to a magnetization component in addition to, and oriented at an angle to the ferromagnetic order in the system. In contrast, an external magnetic field B applied perpendicularly to the ferromagnetic order can induce an electric dipole moment via a mechanism resembling the Coulomb force, where the scalar potential is replaced by the vector potential for B.
Magnetoelectricity occurs most prominently in antiferromagnetically ordered materials, where an electrically induced magnetization is not masked by an intrinsic magnetization in the system. Similar to ferromagnetic order, antiferromagnetic order can have a spin component and an orbital component, and we can have spontaneous antiferromagnetic order due to a staggered exchange field in the material, but the order can also be induced in both paramagnets and ferromagnets. Here we consider the antiferromagnetic diamond structure shown in Fig. 1(c). To study the magnetoelectricity exhibited in quantum wells made from such a material, we develop a k · p envelope-function theory for itinerant-electron diamond antiferromagnets, which is in itself an important result presented in this work. On the basis of this theory, we are able to define an operator τ in terms of itinerantelectron degrees of freedom such that a nonzero expectation value τ signals collinear antiferromagnetic order Zincblende structure that breaks inversion symmetry. (c) Antiferromagnetic diamond structure that breaks time reversal symmetry Θ and inversion symmetry I (though the joint operation ΘI remains a good symmetry). Materials with structure (a) are not magnetoelectric. Those with structure (b) become magnetoelectric when they are magnetized, whereas materials with structure (c) are intrinsically magnetoelectric.
in the same way that a nonzero expectation value σ of the charge carriers' spin operator σ signals ferromagnetic order of spins. Applying our theoretical framework to antiferromagnetically ordered quantum wells placed into external magnetic and electric fields, we reveal them to exhibit magnetoelectric couplings remarkably similar to those found for the ferromagnetically ordered zincblende quantum wells described above. The magnetoelectric response of the antiferromagnetic system can be related to the modification of the quadrupolar equilibrium-current distribution associated with antiferromagnetic order by external electric and magnetic fields.
Analytical results obtained from effective two-band models of confined charge carriers elucidate the basic physical phenomena associated with magnetoelectricity in para-, ferro-and antiferromagnetic quantum wells. Accurate numerical calculations utilizing realistic 8 × 8 and 14×14 k·p Hamiltonians establish a typically large, practically relevant magnitude of the electric-field-induced magnetization in hole-doped quantum wells made from zincblende ferromagnets or diamond-structure antiferromagnets. The ability to illustrate the full complementarity of magnetoelectric responses within the same microscopic theory distinguishes our approach from most previous ones [13]. Our findings provide a platform for further systematic studies aimed at manipulating charges, currents, and magnetic order in solids.
The remainder of this Article is organized as follows. In Sec. II, we define the relevant quantities of interest for our study, establishing the relation between the thermodynamic definitions of polarization (2a) and magnetization (2b) and the electromagnetic definitions of these quantities. We then proceed, in Sec. III, to calculate magnetoelectric responses of quasi-2D electron and hole systems realized in zincblende heterostructures having a Zeeman spin splitting due to an external magnetic field or due to the coupling to ferromagnetic exchange fields. In Sec. IV, we develop a general framework for the k · p envelope-function description of antiferromagnetic order. We use this framework to perform a comprehensive analysis of magnetoelectric phenomena in quantum wells made from diamond-structure antiferromagnets. We summarize our conclusions and provide a brief outlook in Sec. V. Ancillary results are presented in the Appendices.

II. ELECTRIC AND MAGNETIC RESPONSES IN QUASI-2D SYSTEMS
We consider a quasi-2D system in the (x, y) plane with open boundary conditions in the z direction in the presence of a perpendicular electric field E z and an in-plane magnetic field B = (B x , B y ). Throughout this work, vectors like B that have only in-plane components will be indicated by a subscript ' ', and their vanishing z component will be suppressed. Very generally, the polarization and magnetization can be obtained from the free-energy density F via the relations [11] More accurately, the polarization and magnetization only depend on the change of the free energy δF ≡ F (E z , B ) − F (0, 0) due to the fields E z and B .
To simplify the analysis, we assume that only the itinerant charge carriers in the quasi-2D system contribute to the electric and magnetic response. We assume that the confining potential V (z) of the quasi-2D system includes the electrostatic potential due to compensating charges and external gates that ensure overall charge neutrality and that are assumed to be fixed in space. Also, we assume that the potential V (z) defining a quantum well for the quasi-2D system is symmetric, i.e., V (−z) = V (z). We denote the Hamiltonian for the charge carriers by H. The electric field E z enters H via the additional potential eE z z, where e ≡ |e| is the elementary charge. The magnetic field B enters H via the vector potential A that is related to the magnetic field via B = ∇ × A, with ∇ denoting the gradient w.r.t. the position vector r ≡ (x, y, z). In addition, B may enter H via a Zeeman term (g/2)µ B σ · B , where g denotes the g factor, µ B ≡ e /(2m 0 ) is the Bohr magneton, with m 0 being the mass of free electrons, and σ is a dimensionless spin operator [40]. The eigenstates of H associated with eigenvalues E nk have the general form Here n labels the quasi-2D subbands, and k ≡ (k x , k y ) is the in-plane wave vector. The free-energy density can then be written in the form where is the 3D number density, and N s = dz ρ(z) is the 2D (sheet) density of charge carriers in the quantum well. Thus we can rewrite the polarization (8) as where P 0 ≡ −eN s , and the dimensionless number P = z /w describes the average polarization per particle. Similar to the polarization P z , the magnetization M is also the sum of two contributions, with which again represent the energetic and the quantumkinetic effects of B , respectively. Given that B generally enters the Hamiltonian H via both the vector potential A and also via the Zeeman term, the contribution M e can be split further into orbital and spin contributions, where To obtain Eqs. (16a) and (16b), we used once again the Hellmann-Feynman theorem. The first term M o represents the in-plane orbital magnetization [43,44]. In Eq. (16a), the symbol v ≡ ∂H/(∂ k ) denotes the in-plane component of the velocity operator, and with {A, B} ≡ 1 2 (AB + BA). The expression (16a) is associated with the vector potential A = z B ×ẑ that is adopted throughout our work as the appropriate gauge for quasi-2D systems. This is the reason why Eq. (16a) differs from the conventional formula for the orbital magnetization [44] that is obtained for the symmetric gauge [45] A sym ≡ 1 2 B × r, see Appendix A. Similar to P e z , the magnetization M o of a quasi-2D system avoids the technical problems inherent in studies of the bulk (3D) orbital magnetization [44]; it is unambiguously defined independent of the origin of the coordinate system. An orbital magnetization M o is generally accompanied by a nonvanishing in-plane current distribution though in thermal equilibrium, the total current J = dz j (z) is always zero. These currents j (z) are nondissipative because they are not driven by an electric field. (Throughout this work, we assume E = 0 for the in-plane electric field.) Direct experimental observation of the currents j (z) seems impossible, as their nature appears to preclude any ability to make contact to them. However, their ramification in terms of the magnetization M o is detectable.
The second term S in Eq. (15) represents the spin magnetization, given in Eq. (16b) in terms of the dimensionless spin polarization σ nk of individual states. We rewrite this as where S is the dimensionless average spin polarization of the entire system. Similarly, it is convenient to define A polarization P e z represents the dipole term (l = 1) in a multipole expansion of a charge distribution ρ(z) [43]. Similarly, an orbital magnetization M o represents the dipole term (l = 1) in a multipole expansion of a current distribution j (z). Charge neutrality of a localized charge distribution ρ(z) generally requires a vanishing monopole (l = 0) for the multipole expansion of ρ(z). Similarly, a localized current distribution j (z) requires a vanishing monopole for the multipole expansion of j (z). An equilibrium current distribution j (z) that breaks time-reversal symmetry is permitted in ferromagnets and in antiferromagnets [46]. The finite magnetization in ferromagnets implies that the equilibrium current distribution j (z) includes a dipolar component, whereas the vanishing magnetization in antiferromagnets requires equilibrium currents to be (at least) of quadrupolar type. In the context of the present work, we indeed find that the equilibrium current distribution j (z) in antiferromagnetic quantum-well structures gives rise to a quadrupole moment In the quasi-2D geometry relevant for our purpose, the nonzero components of Q defined within classical physics are [47,48] A quantum-mechanical description of the quadrupole moment needs to take into account that, generally, the operators z and v = (v x , v y ) do not commute. At minimum, two definitions for an operator Q are plausible [49], and These two definitions for the operator Q yield different results if the operator v is at least of second order in k z . We are not aware of a satisfactory resolution for this ambiguity. We avoid this problem in the present study by evaluating Q using Hamiltonians H for which v = ∂H/(∂ k ) contains terms that are at most linear in k z . For finite systems, the lowest nonvanishing multipole in a multipole expansion is generally independent of the origin of the coordinate system and in that sense welldefined, whereas higher multipoles depend on the choice for the origin [43]. We therefore limit our discussion below to the lowest nonvanishing multipole. As mentioned above, in infinite periodic crystals, even the lowest nonvanishing multipole moment requires a more careful treatment [44].
As integrals can be more easily and more reliably calculated numerically than derivatives [50], it is more straightforward to evaluate numerically the integrals defining the electromagnetic parts P e z , M o , and S of the response functions. On the other hand, it is more difficult to evaluate accurately the full response functions P z and M that require a numerical differentiation of the free energy F as a function of the applied external fields [51,52]. A detailed account of these technical issues is beyond the scope of the present work. In the following, we thus focus on P e z , M o , and S alone. This is adequate for scenarios where the quantum-kinetic parts P q z and M q of the response functions are less important, which we have found to be generally the case for a strong confinement V (z).

III. MAGNETOELECTRICITY IN ZINCBLENDE PARAMAGNETS AND FERROMAGNETS
A. The model The diamond crystal structure is shown in Fig. 1(a). Space inversion is a good symmetry in diamond so that electronic states are at least twofold degenerate throughout the Brillouin zone [37]. The diamond structure is realized in group-IV semiconductors including C, Si, and Ge. In a zincblende structure, the atomic sites in a diamond structure are alternatingly occupied by two different atoms such as Ga and As or In and Sb [ Fig. 1(b)]. Thus spin degeneracy of the electronic states is lifted in paramagnetic zincblende structures except for k = 0. Spontaneous ferromagnetic order is realized in semiconductors with zincblende structure such as GaMnAs [53] and InMnSb [54], where the ferromagnetic coupling between local Mn moments is mediated by itinerant holes [38,39].
For common semiconductors with a zincblende structure, such as GaAs, InAs, and InSb, the electronic states in a quantum well can be described by a multiband Hamiltonian [55] Here H k is the inversion-symmetric part of H, and H D subsumes Dresselhaus terms due to bulk inversion asymmetry (BIA). V (z) is the quantum-well confinement, so that the wave vector k = (k x , k y ) is a good quantum number, whereas k z becomes the operator −i∂ z . An external electric field E z can be included in H by adding the potential eE z z. Similarly, an external in-plane magnetic field B can be included in H via the vector potential A = z B ×ẑ. In H k + H D we then replace k by the kinetic wave vector k = k + e A. The Zeeman term H Z includes contributions from both the external field B and possibly a ferromagnetic exchange interaction represented by an internal exchange field X that is likewise assumed to be in-plane. A finite exchange field X corresponds to a finite spontaneous magnetization M s in the expansion (1). For X = 0, the system is a paramagnet, where the lowest-order term in the expansion (1) that de- signifying the fact that the system's magnetization scales with the applied field B. For the magnetoelectric effect studied here, a finite Zeeman term H Z indicates, first of all, a breaking of time-reversal symmetry so that the origin of H Z is largely irrelevant for the microscopic mechanism yielding the magnetoelectric response. Nonetheless, as to be expected, we will see below that only for X = 0, the final result for the lowest-order magnetoelectric contribution to the free energy (1) can be expressed via a tensor α ij , whereas in paramagnets the linear dependence of H Z on B is the reason why in lowest order we get terms in Eq.
(1) that are weighted by a third-rank tensor β ijk .
The diagonalization of the Hamiltonian (24) yields the eigenenergies E nk with associated bound states Φ nk (z) ≡ z|nk , where n is the subband index. In the numerical calculations presented below, we use for H the 8 × 8 Kane model and the 14 × 14 extended Kane model as defined in Table C.5 of Ref. [55]. Confinement in the quasi-2D system is due to a finite poten- The numerical solution of H is based on a quadrature method [56]. We evaluate k-space integrals such as Eq. (8) by means of analytic quadratic Brillouin-zone integration [57].
Before presenting numerical results for multi-band models, we illustrate the physical origin and ramifications of magnetoelectricity in zincblende-semiconductor quantum wells by analytical calculations. Specifically, we consider a 2 × 2 model for the Γ 6 conduction band with where m denotes the effective mass, H D is the Dresselhaus term with prefactor d, cp denotes cyclic permutation of the preceding term, σ ≡ (σ x , σ y , σ z ) is the vector of Pauli matrices, and H Z is the Zeeman term that depends on the total field Z ≡ (g/2) µ B B + X . The relation between the simplified Hamiltonian H and the more complete Hamiltonian H is discussed in more detail, e.g., in Ref. [55]. From now on, the direction of Z is chosen as the spinquantization axis for convenience. We will be interested in terms at most quadratic in k and linear in B , where the latter is justified for weak fields B , i.e., when the well width w is smaller than the magnetic length /|eB |. (27) and ϕ Z is the angle between the total Zeeman field Z and the crystallographic direction [100]. The usefulness of writing H as in Eq. (26b) will become clear later on. For E z = 0 and B = 0, the Hamiltonian is

Then the Hamiltonian H becomes [58]
with The eigenstates of with k 2 z = ν|k 2 z |ν . For Z = 0, the spectrum E νσ,k satisfies time reversal symmetry, E νσ,−k = E νσ,k . For Z = 0, the relation E νσ,−k = E νσ,k reflects broken time-reversal symmetry. The latter is a prerequisite for the magnetoelectric effect, as discussed above.
Figures 2(a) and 2(c) illustrate the dispersion (30) for a quasi-2D electron system in a ferromagnetic InSb quantum well with X x = 8 meV, width w = 150Å, and with an electron density N s = 1.0 × 10 11 cm −2 . The numerical calculations in Fig. 2 are based on the more accurate multiband Hamiltonian H introduced above. Band parameters for InSb are taken from Ref. [55].

B. E-induced magnetization
To calculate the equilibrium magnetization induced by an electric field E z , we start from the Hamiltonian (26). Specializing to B = 0 yields with D given by Eqs. (29a), (29b), and (29c), respectively, and Treating the electric field E z in first-order perturbation theory, the eigenstates become with expansion coefficients It will be seen below that, for the calculation of the electric-field-induced magnetization, we can ignore the D that yields an effect of higher order in the Dresselhaus coefficient d. In the following, . . . denotes the average in the unperturbed state |ν , whereas ⟪ . . . ⟫ denotes the average in the perturbed state |νσ (1) in the presence of the electric field E z .
For the equilibrium magnetization (16a), we need to evaluate expectation values ⟪ {z , v (k )} ⟫ using the ve-locity operator associated with the Hamiltonian (31) v We get Up to Eq. (36c), the steps are exact in the sense that they do not assume a perturbative treatment of H M← E . To obtain Eq. (36d), we exploited the fact that the eigenstates |ν of the unperturbed problem can be chosen such that all matrix elements in Eq. (36) become real. For the last line of Eq. (36), we assumed that the potential V (z) is symmetric. The first term in Eq. (36e) yields a vanishing contribution when summed over the equilibrium Fermi sea, as it is proportional to the system's total equilibrium current. Therefore, a nonzero magnetization is due to the second term in Eq. (36e), which yields a contribution independent of the wave vector k . Summing over the Fermi sea and assuming a small density N s such that only the lowest subband ν = 0 is occupied, we obtain for the magnetization (16a) [58] with where l d ≡ 2m 0 d/ 2 is the length scale associated with Dresselhaus spin splitting [59], and distinguishes between a partially and a fully spinpolarized system. For ϕ Z = nπ/2 (n integer), the E zinduced magnetization M o is oriented perpendicular to the field Z. More generally, a clockwise rotation of Z implies a counterclockwise rotation of M o .
The value obtained for the sum in Eq. (38) depends on particularities of the quantum-well confinement. Peculiarly, the sum vanishes for a parabolic (i.e., harmonicoscillator) potential. In contrast, assuming an infinitely deep square well of width w, we get Figure 3(a) illustrates the E z -induced orbital magnetic moment per particle for a ferromagnetic InSb quantum well with width w = 150Å and electron density N s = 1.0 × 10 11 cm −2 . Results in Fig. 3 are based on the more accurate multiband Hamiltonian H. The magnetization (37) complements the more trivial magnetization M tot Z = S Z + M Z that we get already in the absence of a field E z , which is oriented (anti)parallel to Z. The spin magnetization S Z is due to an imbalance between spin eigenstates induced by the Zeeman term (29b). The orbital magnetization M Z is due to spin-orbit coupling. Just like S Z , the orbital contribution is already present in inversion-symmetric diamond structures, i.e., it is a manifestation of spin-orbit coupling beyond the Dresselhaus term (29c) and beyond the simple 2×2 model studied in this section. Therefore, M Z is always present in the numerical calculations based on H. An analytical model for M Z based on H is discussed in Appendix B.
It is illuminating to relate the magnetization (37) to the equilibrium current distribution (18). Using φ ν (z) ≡ z|ν , the perturbed wave functions read In the following, we suppress the argument z of φ ν for the sake of brevity . Using the velocity operator (35), we get in first order of E z and d with matrix elements (spin σ = +) In thermal equilibrium, the first term in Eq. (42b) averages to zero in Eq. (18a). The remaining terms are independent of k so that, for Z = 0, they do not average to zero in Eq. (18a).
The matrix elements contributing to the second term in Eq. (42b) are nonzero independently of an electric field E z (provided the product ν ν is also even). For ν = 2, we get equilibrium currents proportional to φ 0 (z)φ 2 (z) that give rise to a magnetic quadrupole Q [Eq. (21)]. The quadrupolar currents are illustrated in numerical calculations for a quantum well with finite barriers and using the more complete multiband Hamiltonian H, see Figs. 4(a) and 4(c). The quadrupolar currents and the magnetic quadrupole Q are indicative of orbital antiferromagnetic order that is induced parallel to the Zeeman field Z by the interplay of Z, the Dresselhaus term (29c), where we assumed, as before, that only the lowest subband ν = 0 is occupied. As we have τ Z we can interpret such a scenario as ferrimagnetic order. This classification proposed here applies, in particular, to Mndoped semiconductors such as GaMnAs and InMnSb [60]. It is a peculiarity of an infinitely deep potential well that κ ν ν ∝ δ ν ν so that within this model we do not obtain quadrupolar equilibrium currents and orbital antiferromagnetic order.
The last term in Eq. (42b) (with ν = 1) describes E z -induced dipolar currents that contribute to the magnetization [ Fig. 4(c)]. For ϕ Z = nπ/2 (n integer), the quadrupolar and dipolar currents flow (anti)parallel to the field Z, consistent with Eq. (37). As to be expected, the total current J = dz j (z) always vanishes. The fact that the coupling of the currents j (z) to a perpen-dicular electric field E z is dissipationless resembles the Lorentz force. However, it needs to be emphasized that the equilibrium currents j (z) and their manipulation via electric fields are pure quantum effects with no classical analogue.
The numerical calculations for a ferromagnetic quantum well based on the multiband Hamiltonian H and presented in Figs. 4(a) and 4(c) assume that the exchange field X is oriented in x direction. In this case, the equilibrium currents j (z) represented by Eq. (42b) are oriented likewise in x direction. These currents are complemented by equilibrium currents j y representing the orbital magnetization M Z induced by the exchange field X z and discussed in more detail in Appendix B.

C. B-induced electric polarization
To calculate the equilibrium electric polarization induced by a magnetic field B , we start again from the Hamiltonian (26). Specializing to E z = 0 yields with H (0) , H where ϕ B is the angle between the direction of the applied magnetic field B and the [100] crystallographic direction. The perturbation H P←B yields the perturbed states Here, the first term z νσ vanishes for a symmetric potential V (z). The first term in the square brackets describes a k -dependent shift [61][62][63][64] that yields a vanishing contribution to P e z when summed over the equilibrium Fermi sea. Therefore, a nonzero polarization is due to the second term in the square brackets, which yields a contribution independent of the wave vector k . Summing over the Fermi sea, we obtain [58] where λ d is given by Eq. (38). We see that the induced magnetoelectric effects are most pronounced when . In paramagnetic systems with X = 0 and Z = (g/2)µ B B , we have ζ = ±1 when B [110]. In this case the polarization P e z depends quadratically on B , consistent with Eq. (2). Thus the system exhibits a higher-order magnetoelectric effect [23,24] that is the nondissipative counterpart of the previously discussed magnetically induced electric polarization in a multi-quantum-well system [65][66][67]. Figure 3(c) illustrates the polarization (49) for a ferromagnetic InSb quantum well.
The mechanism for the B-induced polarization can be understood as follows: the vector potential A of a magnetic field B has previously been used as a tool to manipulate the charge density ρ(z) in quasi-2D systems such as semiconductor quantum wells. Ordinarily, a field B makes the charge distribution ρ(z) bilayerlike by pushing ρ(z) towards the barriers, but ρ(z) still preserves the mirror symmetry of a symmetric quantum well [61][62][63][64]. This effect stems from terms quadratic in A that we have ignored in the above analytical model. In a low-symmetry configuration [indicated here by the presence of the Dresselhaus term H D given in Eq. (25c)], odd powers of the vector potential A can change ρ(z) in a way that no longer preserves the mirror symmetry of the confining potential V (z). This effect resembles the Coulomb force, where the scalar potential is replaced by the vector potential A. However, it needs to be emphasized that, similar to Landau diamagnetism, we have here a pure quantum effect; it has no classical analogue. This effect is orbital in nature; it does not require a spin degree of freedom. For example, it exists also in spinless 2D hole systems that have a purely orbital Dresselhaus term.
where we ignored terms O(E 2 z ) and O(B 2 ). When averaging over all occupied states, the terms ∝ A drop out. Using Eq. (38), we get [58] consistent with Eqs. (37) and (49). Hence, within the present model, we have P z = P e z and M = M e . The expression (52) can be written as a sum of terms of the type appearing in the third line of the general expansion (1). To illustrate this, we consider the case Z < π 2 N s /m and find Clearly, α ij = 0 requires spontaneous ferromagnetic order due to a finite exchange field X , whose presence will generally also facilitate higher-order terms of the type ∝ γ ijk in Eq. (1). In contrast, β ijk = 0 occurs even in paramagnets, which is consistent with basic symmetry considerations [23,24,[27][28][29] as our zincblendestructure-based system of interest allows for piezoelectricity.
The magnetoelectric contribution (52) to the free energy can also be expressed as in terms of a magnetoelectric vector [68]  The angular dependence of the magnetoelectric effect is governed by the orientation of the vectorτ , which in turn is determined by the orientation of the Zeeman field Z. Indeed, we have a one-to-one correspondence between the orientation of the vector Z in position space and the vector k 0 in reciprocal space. More specifically, the part k (z) 0 of k 0 proportional to σ z turned out to be relevant for the magnetoelectric effect in the above analysis. This vector k (z) 0 is (anti)parallel to the vectorτ . Figure 5(a) shows the relation between the orientation of Z and the orientation of k (z) 0 τ . A similar pattern exists in the Edelstein effect for systems with Dresselhaus spin-orbit coupling (25c) for the orientation of the induced spin polarization as a function of the orientation of an in-plane electric field [55].
The vector E zẑ × B in Eq. (54a) is a toroidal vec-tor, i.e., it is odd under both space inversion and time reversal [70]. On the other hand, Eq. (54b) shows that the vectorτ transforms like a magnetic field, i.e., it is even under space inversion and odd under time reversal. The different transformational properties of the vectors E zẑ × B andτ in Eq. (54a) reflect the broken space inversion symmetry in a zincblende structure. The term δF ∝ E z B in Eq. (52) is generally complemented by a second magnetoelectric term ∝ E z B . This is because the Hamiltonian H also includes a term characterizing the bulk zincblende structure that underlies the quasi-2D systems studied here. The prefactor b is given in Eq. (7.5) of Ref. [55] in terms of momentum matrix elements and energy gaps appearing in the larger Hamiltonian H, yielding b = −221Å/eV for InSb and −1.36Å/eV for GaAs. The term (55) produces a second magnetoelectric term in the free energy, that complements δF in Eq. (52). Their ratio is given by where the expression on the far r.h.s. of Eq. (57) is obtained using Eq. (38) for a hard-wall confinement V (z). This ratio evaluates to 9300/(w[Å]) 2 in InSb and 330/(w[Å]) 2 in GaAs, and it is consequently much smaller than 1 for typical quantum-well widths w 150Å. Experimental signatures of an E z B -type magnetoelectric coupling have recently been observed for charge carriers in deformed donor bound states [71].

E. Magnetoelectricity in ferromagnetic hole systems
The magnetoelectric response obtained in the realistic calculations for electron systems in ferromagnetic InSb quantum wells is small [ Figs. 3(a) and 3(c)]. The response can be greatly enhanced by a suitable engineering of the band structure E nk of the quasi-2D systems. Here quasi-2D hole systems are long known as a versatile playground for bandstructure engineering, where the dispersion of the first heavy-hole (HH) subband is strongly affected by the coupling to the first light-hole (LH) subband [55,72,73]. Figures 6(a) and 6(d) illustrate this for quasi-2D hole systems in paramagnetic InSb quantum wells with width w = 150Å [ Fig. 6(a)] and w = 300Å [ Fig. 6(d)], where HH-LH coupling results in a highly nonparabolic dispersion E 0k of the doubly degenerate ground HH subband. Furthermore, the dispersion is also highly anisotropic, which reflects the cubic symmetry of the underlying crystal structure. An important aspect for the magnetoelectric response is the breaking of time reversal symmetry so that E n,−k = E nk . The interplay between a ferromagnetic exchange field X and HH-LH coupling can result in a highly asymmetric band structure of quasi-2D HH systems with multiple disconnected parts of the Fermi surface, as illustrated in Figs. 6(b) and 6(e) for ferromagnetic InSb quantum wells [54]. Figures 7(a) and 7(c) exemplify the E z -induced orbital magnetic moment per particle, which can rise as high as ∼ 1 µ B for moderate electric fields E z . Figures 8(a) the equilibrium currents. Finally, Figs. 9(a) and 9(c) show the B -induced displacement z that represents the electrostatic polarization via Eq. (12). Unlike the electron case discussed above, the hole systems show a strongly nonlinear dependence of the magnetoelectric response as a function of the applied fields. This result can be understood from the dispersions depicted in Figs. 6(b) and 6(e). With increasing fields, the disconnected parts of the Fermi sea that are located away from k = 0 get depopulated and eventually disappear. The field-induced response drops again when finally only the central part of the Fermi sea around k = 0 accommodates all charge carriers.

IV. MAGNETOELECTRICITY IN DIAMOND ANTIFERROMAGNETS
Space inversion symmetry of a diamond structure is broken in the zincblende structure [ Figs. 1(a) and 1(b)]. Opposite magnetic moments placed alternatingly on the atomic sites of a diamond structure result in an antiferromagnetic structure [ Fig. 1(c)]. Both time reversal Θ and space-inversion symmetry I are broken in such a diamond antiferromagnet. The joint operation ΘI, however, remains a good symmetry so that, similar to paramagnetic diamond, a two-fold spin degeneracy is preserved throughout the Brillouin zone. Nonetheless, as these symmetries are broken individually, invariants proportional to the Néel vector N appear in the Kane Hamiltonian that are forbidden in paramagnetic systems because of time-reversal symmetry. These invariants are derived in Sec. IV A. The diamond structure is realized by the A atoms of spinels AB 2 X 4 . Frequently, spinels with magnetic A atoms give rise to highly frustrated magnetic order [74]. Beyond that, a recent study combining experiment and theory [75] identified CoRh 2 O 4 as a canonical diamondstructure antiferromagnet, where the magnetic moments on the Co atoms are ordered as shown in Fig. 1(c) (magnetic space group I4 1 /a m d).

A. The model
Our goal is to incorporate the effect of antiferromagnetic order into the k·p envelope-function theory [55,76] underlying multiband Hamiltonians as in Eq. (24). To this end, we start from the well-known sp 3 tight-binding model for diamond and zincblende structures with spinorbit coupling included [77,78]. This model includes the s-bonding valence band Γ v 6 , the p-bonding valence bands Γ v 8 and Γ v 7 , the s-antibonding conduction band Γ c 6 , and the p-antibonding conduction bands Γ c 8 and Γ c 7 . Except for the low-lying valence band Γ v 6 , these bands are also the basis states for the 14 × 14 extended Kane model [55,79].
We add a staggered exchange field Y on the two sublattices of the diamond structure as depicted in Fig. 1(c). Using the phase conventions for these basis functions of H that are described in detail in Appendix C of Ref. [55], the field Y yields terms in the off-diagonal blocks H 8c 8v , H 8c 7v , and H 7c 7v of H that are listed in Table I using  the notation of Table C.5 in Ref. [55]. The vector N denotes the Néel unit vector with components  Fig. 1(c) appear already for k = 0. In the diagonal blocks H 6c 6c , H 8v 8v , and H 7v 7v , a Taylor expansion of the tightbinding Hamiltonian about k = 0 yields mixed terms proportional to powers of components of Y and powers of components of k. The lowest-order invariants obtained in this way are also listed in Table I. Alternatively, these terms can be derived by means of quasi-degenerate per-  [55], and N i ≡ Yi/Y denotes the Cartesian components of the unit vector parallel to the staggered exchange field Y on the sublattices of the diamond structure [see Fig. 1(c)].
turbation theory [55] applied to H with H Y 8c 8v , H Y 8c 7v , and H Y 7c 7v included. The latter approach yields explicit, albeit lengthy, expressions for the prefactors d and D i jj as a function of Y = | Y| that are omitted here. As to be expected for antiferromagnetic diamond, the Ydependent invariants in Table I break time-reversal symmetry, but they do not lift the spin degeneracy. Using quasi-degenerate perturbation theory, we also obtain several invariants in the valence band block H Y 8v 8v that are proportional to both Y and an external electric field E. These invariants are listed in Table I as well. They describe a spin splitting proportional to the field E (but independent of the wave vector k) that is induced by the antiferromagnetic exchange coupling. All invariants listed in Table I can also be derived by means of the theory of invariants [76] using the fact that the staggered exchange field Y is a polar vector that is odd under time reversal.
According to Table I, in lowest order the Γ 6 conduction band in a diamond antiferromagnet is described by the Hamiltonian with H k given in Eq. (25b), and where d is a prefactor proportional to Y. Formally, H N has the same structure as the Dresselhaus term (25c), with the spin operators σ i replaced by the numbers N i and d replaced by d. Therefore, the following study of magnetoelectric coupling in antiferromagnetic diamond proceeds in remarkable analogy to the study of magnetoelectric coupling in a paramagnetic or ferromagnetic zincblende structure presented in Sec. III [69]. As H N and, in fact, the entire Hamiltonian (58a), do not de-pend on the the charge carriers' spin, the latter will be a silent degree of freedom in the following considerations. For the analytical model studied below, it is easy to see that a purely in-plane Néel unit vector N yields the largest magnetoelectric coupling. Assuming therefore that N is oriented in-plane, the full Hamiltonian becomes [including terms up to second order in k , compare Eq. (26)] Here ϕ N denotes the angle that N makes with the x axis, and we introduced the operator For E z = B = 0 and treating H N in first order, the subband dispersions become which are parabolae that are shifted in the (k x , k y ) plane by k 0 ν . The shift k 0 ν is a fingerprint for the broken time-reversal symmetry in the antiferromagnet.  Table I. One can approximate these results with the smaller Hamiltonian (59) using d ≈ 80 eVÅ 3 .

B. The Néel operator
We now digress to discuss a few general properties of the model for antiferromagnetic order proposed here. It is well-known that the Zeeman term (25d) with an exchange field X provides a simple mean-field model for itinerant-electron ferromagnetism. Similarly, H N is a phenomenological model for collinear itinerant-electron antiferromagnetism.
The operator conjugate to the ferromagnetic exchange field X is the (dimensionless) spin-polarization operator σ = ∂H/∂X. In the mean-field theory underlying the present work, a nonzero expectation value σ indicates ferromagnetic order of spins. Similarly, the operator conjugate to the staggered exchange field Y is the (again dimensionless) Néel operator for the staggered magnetization, where the prefactor q τ ≡ d/Y depends on the momentum matrix elements and energy gaps characterizing the Hamiltonian H, but it is independent of the exchange field Y. A nonzero expectation value τ indicates collinear orbital antiferromagnetic order. Like the staggered exchange field Y, the Néel operator τ is a polar vector that is odd under time reversal. Thus τ = 0 represents a (polar) toroidal moment [68,70,80]. On the other hand, X and σ are axial vectors that are odd under time reversal. In that sense, σ and τ quantify complementary aspects of itinerant-electron collinear magnetic order in solids.
In systems with spin-orbit coupling such as the ones studied here, the spin magnetization S associated with the expectation value σ is augmented by an orbitalmagnetization contribution, yielding the total magnetization M. A magnetization M arises due to the presence of an exchange field X or an external magnetic field B, but it may also arise due to, e.g., an electric field E (the magnetoelectric effect studied here) or a strain field (piezomagnetism [11,81,82]). Similarly, a nonzero expectation value τ can be due to a staggered exchange field Y. But it may also arise due to, e.g., the interplay of an exchange field X, spin-orbit coupling, and confinement [Eq. (42b)]. The theory for M and τ induced by external perturbations can be phrased very generally using the theory of material tensors [76,83,84]. Each of the induced effects comes commonly in a dissipative and a nondissipative version, similar to the generation of magnetic-moment densities in the dissipative Edelstein effect in paramagnetic systems [85][86][87][88] and the nondissipative magnetoelectric effect in magnetically ordered systems [27][28][29].

C. E-induced magnetization
To calculate the equilibrium magnetization, we start from the Hamiltonian (59). Treating the electric field E z in first-order perturbation theory, the eigenstates become compare Eq. (33). For the equilibrium magnetization (16a), we need to evaluate expectation values ⟪ {z , v (k )} ⟫ using the velocity operator associated with the Hamiltonian (59) We get compare Eq. (36). Again, we ignored any k or k 0 dependence of the perturbed states |ν (1) , which is a higherorder effect. The first term in Eq. (65e) yields a vanishing contribution when summed over the equilibrium Fermi sea, as it is proportional to the system's total equilibrium current. Therefore, a nonzero magnetization is due to the second term in Eq. (65e), which is independent of the wave vector k . We can obtain Eq. (65e) from Eq. (36e) by replacing ϕ Z with ϕ N and putting σ = 1 for all states. The latter implies that the effect described by Eq. (65e) is maximized compared with Eq. (36e) because both spin orientations in the antiferromagnet contribute constructively. Summing over the Fermi sea, we obtain for the magnetization (16a) with and l d ≡ 2m 0 d/ 2 , in complete analogy with Eqs. (37) and (38) [69]. For ϕ N = nπ/2 (n integer) the induced magnetization is oriented perpendicular to the Néel vector N. More generally, a clockwise rotation of N implies a counterclockwise rotation of M o . Figure 3(b) illustrates the E z -induced magnetization for an antiferromagnetic InSb quantum well with width w = 150Å and electron density N s = 1.0 × 10 11 cm −2 .
Again, it is illuminating to compare Eq. (66) with the equilibrium current distribution (18). Using φ ν (z) ≡ z|ν , the perturbed wave functions read Using the velocity operator (64), we get in first order of E z and d where the matrix elements κ ν ν are given by Eq. (43) with ϕ Z replaced by ϕ N . Equation (69) is obtained from Eq. (42b) by putting σ = +1 so that the interpretation of Eq. (69) proceeds similarly. In thermal equilibrium, the first term in Eq. (69) averages to zero in Eq. (18a). The remaining terms are independent of k so that they do not average to zero in Eq. (18a). The second term (ν = 2) describes a quadrupolar equilibrium current proportional to φ 0 (z)φ 2 (z) independent of the electric field E z . Such quadrupolar orbital currents are a generic feature of antiferromagnets; they are the counterpart of dipolar orbital currents representing the orbital magnetization in ferromagnets (see Appendix B) [89]. Similar to Eq. (44), the orbital antiferromagnetic order can be quantified using the Néel operator τ . The Hamiltonian (59) (with E z = 0) yields As to be expected, we have τ N. The last term in Eq. (69) (ν = 1) describes E zinduced dipolar currents, i.e., a magnetization. In a quantum well of width w, the equilibrium currents j (z, νk ) occur on a length scale of order w, which is typically much larger than the lattice constant of the underlying crystal structure. The magnetic multipoles associated with the current distribution may thus be accessible experimentally. They may even open up new avenues to manipulate the magnetic order in antiferromagnets. Figures 4(b) and 4(d) illustrate the equilibrium currents for antiferromagnetic InSb quantum wells.
It is illuminating to study a second mechanism for an E-induced magnetization based on the antiferromagnetic exchange term (58b) that manifests itself as a spin magnetization (16b). Generally, an electric field E z applied to a quantum well gives rise to a Rashba term [55,90] with Rashba coefficient a R E z , resulting in spin-split eigenstates where ϕ k is the angle between k and the x axis, and we assumed as before that the orbital part |ν of the eigenstates is independent of k . Thus we have Also, Rashba spin-orbit coupling gives rise to an imbalance between the two spin subbands ±, which can be characterized by Fermi wave vectors k F± ≈ √ 2πN s ∓ a R E z m/ 2 . Performing the average (16b) over all occupied states in these spin subbands [assuming a dispersion (61) with small k 0 = 0 and slightly different Fermi wave vectors k F± ], we obtain a nonzero equilibrium spin polarization Inserting this result into (20) yields a spin magnetization that complements the orbital magnetization (66). As to be expected, both terms have the same dependence on the direction of the vector N.
We can interpret the spin polarization (74) as follows. The Rashba term (71) yields a spin orientation (73) of individual states |νk ± . Nonetheless, for nonmagnetic systems in thermal equilibrium, the net spin polarization is zero because time-reversal symmetry implies that we have equal probabilities for the occupation of timereversed states |νk ± and |ν, −k ± with opposite spin orientations. This argument for nonmagnetic systems is closely related to the fact that thermal equilibrium in a time-reversal-symmetric system requires that the Fermi sea is centered symmetrically aboutk = 0. A nonzero shiftk of the Fermi sea, and thus a nonzero average spin polarization, are permitted in nonmagnetic systems as a quasistationary nonequilibrium configuration in the presence of a driving electric field E , which is known as the Edelstein effect [85][86][87][88]91]. The spin polarization (74), on the other hand, is entirely an equilibrium effect. It can occur in antiferromagnetic systems where time-reversal symmetry is already broken in thermal equilibrium as expressed by the shiftk = k 0 .
It follows from Table I that we generally get a spin splitting proportional to E z even at k = 0, which yields a third, Zeeman-like contribution to the total magnetization (20). For quasi-2D hole systems, this effect can be substantial. For quasi-2D electron systems, this effect is of second order in the staggered exchange field Y.

D. B-induced electric polarization
Our goal is to evaluate the polarization (8) The perturbation H P←B yields perturbed states |ν (1) . We get As before [Eq. (48)], the first term z ν vanishes for a symmetric potential V (z). The first term in the square brackets describes a k -dependent shift [61][62][63][64] that yields a vanishing contribution to P e z when summed over the equilibrium Fermi sea. Therefore, a nonzero polarization is due to the second term in the square brackets, which is independent of the wave vector k . Summing over the Fermi sea, we obtain compare Eq. (49) [69]. Figure 3(b) illustrates the Binduced polarization for an antiferromagnetic InSb quantum well.
where we ignored terms O(E 2 z ) and O(B 2 ). When averaging over all occupied states, the terms ∝ A drop out. Using Eq. (67), we get consistent with Eqs. (66) and (77). Decomposing δF into terms present in the third line of Eq. (1) yields with Thus similar to the ferromagnetic case [Eqs. (53)], antiferromagnetic order gives rise to α ij = 0 and could also generate higher-order magnetoelectric contributions of the type ∝ β ijk and ∝ γ ijk in Eq. (1). However, unlike the paramagnetic zincblende structure where β ijk = 0, the high symmetry of a paramagnetic diamond structure precludes the existence of any magnetoelectric effects. Equation (80a) can also be expressed as δF = −τ · (E zẑ × B ) in terms of the magnetoelectric vector which is analogous to the form of the magnetoelectric vector (54b) found for the ferromagnetic case. We havẽ τ k 0 , and, like N, the vectorτ is a toroidal vector. Once again, the nonlinear dependence of the magnetoelectric response on the applied fields is due to the depopulation of the disconnected parts of the Fermi sea that are located away from k = 0.

V. CONCLUSIONS AND OUTLOOK
We present a detailed theoretical study of how magnetoelectricity arises in magnetically ordered quantum wells with broken time-reversal symmetry and broken space-inversion symmetry.
Quasi-2D systems based on zincblende ferromagnets [ Fig. 1(b)] and diamondstructure antiferromagnets [ Fig. 1(c)] exhibit an analogous linear magnetoelectric response, i.e., an in-plane magnetization induced by a perpendicular electric field [Eqs. (37) and (66)], as well as a perpendicular electric polarization arising from an in-plane magnetic field [Eqs. (49) and (77)]. In realistic calculations, the magnitude of the magnetoelectric response is small in quasi-2D electron system (Fig. 3), but it is sizable for quasi-2D hole systems (Figs. 7 and 9). Our findings suggest that bandstructure engineering and nanostructuring are fruitful avenues for generating and tailoring magnetoelectricity in a host of materials.
Both space-inversion symmetry and time-reversal symmetry must be broken so that an electric field can couple to a magnetization and a magnetic field can couple to a polarization. Time-reversal symmetry can be broken in thermal equilibrium due to a spontaneous magnetic order or due to an external magnetic field. Timereversal symmetry is also broken in dissipative nonequilibrium processes. Hence, very generally, two distinct avenues exist for magnetoelectric couplings and related effects. These effects exist in thermal equilibrium in systems with magnetic order. The magnetoelectric effect studied in the present work belongs to this category. On the other hand, when an electric field gives rise to a dissipative current in paramagnetic systems, this can also induce a magnetization, which is known as the Edelstein effect [85][86][87][88] or the kinetic magnetoelectric effect [91]. The detailed requirements for the equilibrium magnetoelectric effect and the nonequilibrium Edelstein effect are very different. The latter exists for those paramagnetic crystal classes that permit axial second-rank tensors [88]. Therefore, the Edelstein effect does not exist, e.g., for the bulk zincblende structure [ Fig. 1(b)], despite the fact that inversion symmetry is broken in the zincblende structure.
Our study yields a new unified picture of magnetic order. Ferromagnetic order is characterized by a magneticmoment density M (a magnetization), which is even under space inversion and odd under time reversal. In itinerant-electron systems, orbital ferromagnetic order is associated with dipolar equilibrium currents. On the other hand, collinear orbital antiferromagnetic order is characterized by a toroidal-moment density τ for the Néel operator τ , which is odd both under space inversion and under time reversal. The toroidal antiferromagnetic order is associated with quadrupolar equilibrium currents. For the itinerant-electron quantum systems studied in the present work, the equilibrium current distributions are slowly varying on the length scale of the lattice constant [Figs. 4 and 8]. The magnetization M and the toroidal-moment density τ quantify complementary aspects of itinerant-electron collinear magnetic order in solids. Ferrimagnetic systems are characterized by both expectation values M and τ being finite simultaneously.
Ferromagnetic order M arises due to the presence of an exchange field or an external magnetic field, but it may also arise due to, e.g., an electric field (the magnetoelectric effect studied here) or a strain field (piezomagnetism [11,81,82]). Similarly, antiferromagnetic order τ can be due to a staggered exchange field. But it may also arise due to, e.g., the interplay of ferromagnetic order, spin-orbit coupling, and confinement [Eq. (42b)]. The theory for how M and τ are induced by external perturbations can be phrased very generally using the theory of material tensors [76,83,84], where each of the induced effects comes generally in a dissipative and a nondissipative version. This view applies, e.g., to recent efforts geared towards an electric manipulation of antiferromagnetic order [92][93][94]: if the symmetry of an antiferromagnetic system is lowered by means of an electric field E, this may result in new nonzero components of the toroidal-moment density τ that are forbidden by symmetry when E = 0. Generally, the manipulation of itinerant-electron ferromagnetic or antiferromagnetic order via external perturbations can be viewed as manipulating the underlying equilibrium current distribution.
To obtain realistic quantitative expressions for magnetoelectricity in quasi-2D systems, we developed a k · p envelope-function theory for charge carriers in antiferromagnetically ordered materials. The explicit form of the Néel operator τ depends on the symmetry of the system under investigation. In the present work, we used the envelope-function theory to derive explicit expressions for τ in antiferromagnetic diamond structures. The theoretical formalism and fundamental understanding of antiferromagnetic order presented in this work can be applied to undertake more comprehensive studies of itinerantelectron antiferromagnets. One area of great relevance [94] is the exploration of antiferromagnetic versions for the Edelstein effect; i.e., possibilities to induce and manipulate antiferromagnetic order using charge currents. More generally, reliable modeling of transport properties for antiferromagnetic-spintronics devices [4,5] requires the level of detail and realism provided by our envelopefunction theory. Basic questions concerning magnetization dynamics in metallic antiferromagnets that are attracting current interest [95] can also be addressed.

ACKNOWLEDGMENTS
RW and UZ acknowledge stimulating discussions with A. Hoffmann and H. Saglam. In addition, they thank A. Hoffmann for support. RW also benefitted from discussions with D. Cahill, M. Gilbert, K. Kang, A. Schleife, M. Shayegan, D. Shoemaker, and J. Sipe. UZ's interest in the magnetoelectricity of quantum wells was initiated by interesting conversations with B. Weber, and he also thanks A. Kamra for useful discussions. This work was supported by the NSF under Grant No. DMR-1310199 and by the Marsden Fund Council (Contract No. VUW1713) from New-Zealand government funding managed by the Royal Society Te Apārangi. Work at Argonne was supported by DOE BES under Contract No. DE-In a single-particle picture for itinerant electrons with kinetic momentum k = k + eA, where A is the vector potential for the magnetic field B = ∇ × A, we get where v = ∂H/(∂ k) = ∂H/(∂ k) is the velocity operator and we took the symmetrized product of noncommuting operators. Repeated indices are summed over. For the symmetric gauge A sym = 1 2 B × r, we have ∂A sym where ijk denotes the totally antisymmetric tensor. Thus which is the conventional formula for the magnetization [44,45] consistent with classical electromagnetism [43].
On the other hand, we get for the asymmetric gauge A = z B ×ẑ employed in the present work whose components differ by a factor of 2 from corresponding terms with r j = z in Eq. (A4). Both expressions for m are consistent with [43] j = −e v = −∇ × m .
Similar to the definition (A1) of the magnetic-moment operator, the operator of the electric dipole moment can be defined as p = −∂H/∂ E. The electric field E can be introduced into H via a scalar potential as in Eq. (32), or via a time-dependent vector potential. Therefore, the explicit form of the electric-dipole moment operator p is also gauge-dependent.

Appendix B: Orbital magnetization induced by Zeeman coupling
Very generally, even in the absence of an electric field E z , the Zeeman term induces a spin magnetization S Z (anti)parallel to the Zeeman field Z and proportional to the g factor. However, in a more complete multiband description similar to the Dirac equation, the g factor for an explicit Zeeman term may be greatly reduced or completely absent. In such an approach, we obtain instead an orbital magnetization M Z due to equilibrium spinpolarized currents, for which spin-orbit coupling plays an essential role. Here we demonstrate that the spin magnetization S Z in a reduced model with g-factor g is equal to the orbital magnetization M Z in the corresponding multiband model.
We illustrate this point for the 8 × 8 Kane model [55,96,97], where the orbital magnetization M Z is due to the off-diagonal coupling between the conduction and valence bands linear in k. The physics that is essential for M Z is thus contained in the simplified Kane Hamiltonian [55] Here E c denotes the conduction band edge (Γ c 6 ), E v ≡ E c − E 0 is the valence band edge (Γ v 8 ) with fundamental gap E 0 , ∆ 0 is the spin-orbit gap between the topmost valence band Γ v 8 and the spin split-off valence band Γ v 7 , and P denotes Kane's momentum matrix element. The terms h c = µ c k 2 z + V c (z) and h v = −µ v k 2 z − V v (z) embody remote-band contributions quadratic in k z with µ c , µ v > 0 and confining potentials V c (z), V v (z) ≥ 0.
While g = 0 for the HamiltonianH, a spin magnetization S Z is obtained whenH is projected on the Γ c 6 conduction band, yielding a 2 × 2 Hamiltonian as in Eq. (25) including a Zeeman term H Z with g-factor g. To express g in terms of the parameters ofH, we decompose H =H (0) +H (1) , whereH (0) contains the diagonal elements ofH, whileH (1) contains the off-diagonal terms linear in k. The eigenstates ofH (0) are bound states |β, νσ ≡ |β, ν ⊗ |σ in the conduction band Γ c 6 (β = c), in the light-hole valence band Γ v 8 (β = l) and in the spin split-off valence band Γ v 7 (β = s) with eigenenergies E β νσ ≡ E β ν + σZ. As before, we introduce an in-plane magnetic field B via the vector potential A = z B ×ẑ. Second-order quasi-degenerate perturbation theory for B then yields Roth's formula [55,98], This calculation is similar to how the g-factor in the Zeeman term of the Pauli equation is derived from the Dirac equation. An imbalance between spin-up and spin-down states (due to an exchange field X or due to an external field B ) thus implies a spin magnetization (16b) proportional to g. For comparison, we now evaluate the orbital magnetization (16a) fromH without projecting on the subspace Γ c 6 . In the following discussion,Z stands for an exchange field X or a magnetic field B that entersH via the vector potential A. Focusing on the states in the conduction band and treatingH (1) in first order perturbation theory, the perturbed eigenstates read |c, νσ (1) = |c, νσ + P where we neglected contributions linear in k as these lead to higher-order corrections in Eq. (B6) below. In the absence of a fieldZ, the eigenstates (B3) are twofold degenerate (σ = ±). The states (B3) are also the appropriate unperturbed states for first-order degenerate perturbation theory for a fieldZ oriented in z direction. If instead, we consider a fieldZ oriented in-plane, the appropriate unperturbed states become |c, νσ, ϕ Z (1) = 1 √ 2 |c, ν+ (1) + σ exp(iϕ Z ) |c, ν− (1) , (B4) where ϕ Z is defined, as before, as the angle betweenZ and the crystallographic direction [100].
The velocity operator is independent of k and independent ofZ. Using the states (B4), the matrix elements (17) of the orbital magnetization can be expressed in the form with g given in Eq. (B2). The matrix elements of the orbital magnetization within the multiband Hamiltoniañ H are thus equal to the matrix elements of the spin magnetization in the two-band Hamiltonian H. In lowest order ofZ, these Hamiltonians yield the same imbalance between the occupation numbers for the respective spin states σ = ±. Thus it follows from Eq. (20b) that, averaged over all occupied states, the orbital magnetization withinH equals the spin magnetization within H. In both approaches, the magnetization vanishes in the limit ∆ 0 → 0. IfZ represents an external magnetic field B , we also get an entirely orbital (spin-independent) diamagnetic contribution to the total orbital magnetization M o [99], which often significantly reduces M o . This compensating diamagnetic shift does not exist in ferromagnets, wherẽ Z represents an internal exchange field X .