Engineering Chiral and Topological Orbital Magnetism of Domain Walls and Skyrmions

Electrons which are slowly moving through chiral magnetic textures can effectively be described as if they where influenced by electromagnetic fields emerging from the real-space topology. This adiabatic viewpoint has been very successful in predicting physical properties of chiral magnets. Here, based on a rigorous quantum-mechanical approach, we unravel the emergence of chiral and topological orbital magnetism in one- and two-dimensional spin systems. We uncover that the quantized orbital magnetism in the adiabatic limit can be understood as a Landau-Peierls response to the emergent magnetic field. Our central result is that the spin-orbit interaction in interfacial skyrmions and domain walls can be used to tune the orbital magnetism over orders of magnitude by merging the real-space topology with the topology in reciprocal space. Our findings point out the route to experimental engineering of orbital properties of chiral spin systems, thereby paving the way to the field of chiral orbitronics.

The field of magnetism is witnessing a recent spark of interest in Berry phase and transport effects which originate in non-collinear magnetism and spin chirality [1][2][3] . One of the recent outstanding observations made in this field is the generation of large current-induced Hall effects in strongly frustrated metallic antiferromagnets 4 and the topological Hall effect (THE) in skyrmions 2,5 . On the other hand, the physics of the fundamental phenomenon of orbital magnetism has been experiencing a true revival which can be attributed to the advent of Berry phase concepts in condensed matter 6,7 . The Berry phase origin of the orbital magnetization (OM) and its close relation to the Hall effect makes us believe that noncollinear spin systems can reveal a rich landscape of orbital magnetism relying on spin chirality rather than spin-orbit interaction (SOI) [8][9][10][11] . The corresponding phenomenon of topological orbital magnetization (TOM) [8][9][10][11] is rooted in the same physical mechanism that drives the emergence of non-trivial transport properties such as the THE in chiral skyrmions or the anomalous Hall effect in chiral antiferromagnets [12][13][14] .
The promises of topological contribution to the orbital magnetization are seemingly very high, since it offers new prospects in influencing and detecting the chirality of the underlying spin texture by addressing the orbital degree of freedom, which is the central paradigm in the advancing field of orbitronics 15 . And while the emergence of topological orbital magnetism in several nm-scale chiral systems has been shown from first principles and tight-binding calculations 8,9 , our understanding of this novel phenomenon is basically absent. In particular, this concerns its conceptually clear definition as well as our ability to tailor the properties of this effect in complex interfacial chiral systems, which often exhibit strong spin-orbit interaction. These are the two central questions we address in this work.
As has been shown in the case of skyrmions, the variety of topological phenomena which arise intrinsically from the non-trivial magnetization configurationn(x, y) can be attributed to an "emergent" magnetic field B z eff 1 . The occurrence of this field is connected to the gauge-invariant Berry phase the electron's wavefunction acquires when traversing the texture [16][17][18] (see Fig. 1b) for an intuitive illustration). In the adiabatic limit, this phase can be attributed to the effect of B z eff , explicitly given by the expression where the sign depends on the spin of the electron. When integrated over an isolated skyrmion, the total flux of B z eff is quantized to integer multiples of 2Φ 0 , where Φ 0 ≈ 2 × 10 3 T nm 2 is the magnetic flux quantum, while the integer prefactor can be identified with the topological charge of a skyrmion, N sk , essentially counting the number of times the spin evolves around the unit sphere when traced along a path enclosing the skyrmion center. Formally, the non-collinear systemn(x, y) can therefore be portrayed as a collinear one, albeit at the price of introducing the magnetic field B z eff into the Schrödinger equation. Just as an ordinary magnetic field would, the emergent magnetic field in chiral systems couples directly to the orbital degree of freedom and provides an intuitive mechanism for the topological Hall effect of skyrmions 19 as well as a possible explanation for the emergence of TOM. Figure 1. Schematic depiction of emergent magnetic fields in a) Néel spirals and b) Néel skyrmions. As electrons (grey spheres) are adiabatically traversing these non-collinear magnetic textures (small arrows, with color indicating the z-projection), their wave function twists just in the same way as it would under the influence of an external magnetic field (the direction is depicted with vertical arrows, the sign and magnitude is illustrated by the colored background). The integrated flux of this emergent topological field over the skyrmion is quantized, while the averaged value of the emergent chiral field for a uniform spin-spiral is zero (although it can be non-zero for a 90 • domain wall). The emergent field locally gives rise to persistent currents (depicted with circular arrows) and the corresponding a) chiral (for a spiral) and b) topological (for a skyrmion) orbital magnetization.
Here, we uncover the emergence of distinct contributions to the orbital magnetization in slowly-varying chiral textures by following the intuition that such contributions should acquire the natural form where χ oms is the orbital magnetic susceptibility of the electronic system 20,21 . Indeed, we demonstrate that in the limit of vanishing SOI the topological orbital magnetization can be expressed in this way. We also discover that in the limit of small, yet non-zero SOI there is a novel chiral contribution to the orbital magnetization described by (2) with the properly defined chiral emergent field which can be finite already for one-dimensional systems (see Fig. 1a)). Moreover, by exploiting a rigorous semiclassical framework, we demonstrate that in interfacial chiral systems with finite SOI, the orbital magnetism can be tuned over orders of magnitude by varying the SOI strength within the range of experimentally observed values. We also underpin the crucial role that the topology of the local electronic structure of textures has in shaping the properties of orbital magnetism in chiral magnets. We discuss the bright avenues that our findings open, paving the way to the experimental observation of this phenomenon and to the exploitation of the orbital degree of freedom in chiral systems for the purposes of chiral orbitronics.

Results
The semiclassical formalism we are referring to in our work is based on the Green's function perturbation theory as presented by Onoda et al. 22 . We put the orbital magnetism of chiral systems on a firm quantum-mechanical ground, formulating a rigorous theory for the emergence of orbital magnetism in non-collinear systems. The motivation for this approach is twofold. First of all, the expression for B eff arises from the adiabatic limit 19,23 , a regime where semiclassical approaches have been successfully applied in order to investigate Berry phase physics 7 . Secondly, this certain type of gradient expansion 24 provides a systematic guide through higher orders of perturbation theory where standard methods would be cumbersome.
It is based on an approximation to the single-particle Green's function and allows us to trace the orders of perturbation theory for chiral magnetic textures, distinguishing corrections to the out-of-plane orbital magnetization 25 M om = 1 M (n)e z of a locally ferromagnetic system which appear as powers of the derivatives of the magnetization with respect to real-space coordinates: where ∂ i = ∂/∂x i . Here and in the following discussion, summation over repeated indices is implied with greek indices α, β ∈ {x, y, z} and latin indices i, j ∈ {x, y}.
The assignment of M tom to the second order expansion, Eq. (4), is based on our intuitive expectation, Eq. (2). The question whether or not this term is "topological" will be discussed in the following and is answered by the semiclassical perturbation theory (see Methods). In addition to the effect of TOM we propose a novel contribution to the orbital magnetization, which is linear in the derivatives of the underlying texture, Eq. (3), and thereby generally sensitive to its chirality. We thus refer to it as the chiral orbital magnetization (COM). We will show how this effect can be attributed to a different kind of effective field (see Fig. 1a)) which emerges from the interplay of spin-orbit coupling and non-collinearity along one spatial dimension.
While our approach is very general, for the purposes of including into consideration the effect of interfacial spin-orbit coupling and providing realistic numerical estimates, we focus our further analysis on the two-dimensional magnetic Rashba model where m * eff is the electron's (effective) mass, σ denotes the vector of Pauli matrices, α R is the Rashba spin-orbit coupling constant, and ∆ xc is the strength of the local exchange field. This model has been proven to be extremely fruitful in unravelling various phenomena in surface magnetism 26 and is known for its pronounced orbital response 27 .
Emergent Fields of Spin Textures. Before discussing the emergence of orbital magnetism in this model, it is rewarding to discuss the appearance of effective fields in slowlyvarying chiral spin textures in the limit of |α R | |∆ xc |. In this regime, it can be shown that to linear order in α R , the spin-orbit coupling can be absorbed into a perturbative correction of the canonical momentum p → p + eA R , with A R ≡ m * eff α R ijz σ i e j /e. This means that the Hamiltonian can be rewritten as: For |α R | |∆ xc | 1 and ∆ xc → ∞ the spin polarization of the wavefunctions is only weakly altered away fromn and we can use an SU (2) gauge field, defined by U † (σ ·n) U ≡ σ z , to rotate our Hamiltonian into the local axis specified byn (neglecting the terms of the order O(α 2 R )) 28,29 : where the potential A now comprises the mixing of two gauge fields: A = U † A R U + A xc , with the additional contribution 1 To be precise, with correct physical dimensions, one should compare the length scales λ R = /α R m * eff and λxc = / ∆xcm * eff and write λ R λxc instead. A xc = −i U † ∇U/e. The essential idea is now the following: as ∆ xc → ∞, electrons are confined to the bands which correspond either to spin-up states |↑ or spin-down states |↓ depending on sgn(∆ xc ). This means that we can effectively replace the vector potential by its adiabatic counterpart, i.e., where A R ad = (U † A R U) ad . Thus, the effective Hamiltonian for ∆ xc → ∞ contains the vector potential of a classical magnetic field which couples only to the orbital degree, accompanying the "ferromagnetic" system 23,30,31 . It is given by the classical expression B eff = ∇ × A ad = B R eff + B xc eff . By following this procedure, one finds the expressions We thus arrive at the fundamental result that in addition to the field given by Eq. (9) above, which can be recognized as the generalization of Eq. (1), there is a contribution to the overall field which explicitly depends on the chirality of the underlying texture and is non-vanishing already for one-dimensional spin textures. In this context, it makes sense to refer to these co-existing fields as topological and chiral for B xc eff and B R eff , respectively, see Fig. (1). Importantly, in contrast to the emergent topological field, (B R eff ) z , the local magnitude of (B R eff ) z is directly proportional to the strength of the spin-orbit interaction as given by α R . This appears to be very promising with respect to achieving a large magnitude of the chiral field in chiral spin textures emerging at surfaces and interfaces. To give a rough estimate, assuming a pitch of the texture on a length scale of L = 20 nm and α R = 1 eVÅ the amplitude of the local chiral emergent field reaches as much as 2πm e α R /(eL) ≈ 270 T, which is roughly by an order of magnitude larger than the corresponding topological field in a skyrmion of a similar size 1 .
The emergence of two types of fields in spin textures, appearing in Eqs. 9 and 10, is crucial for a qualitative understanding of the emergence of topological and chiral orbital magnetism, which are discussed in detail below.

Chiral Orbital Magnetization.
To get a first insight into the novel effect of COM, we consider the limit of small SOI, i.e., α R ∆ xc . In this case, the gradient expansion (see Methods) provides an analytic expression for the local space-dependent orbital moment. Up to O(α R ) it is given by where the function h(x) ≡ (3x 2 − 1)Θ(1 − |x|)/2 describes the energy dependence of COM with Θ representing the Heaviside step-function. The magnitude of COM is thus directly proportional to the strength of spin-orbit interaction and vanishes in the limit of zero α R . Furthermore, M com is proportional to the diamagnetic Landau-Peierls susceptibility 32 χ ↑+↓ LP = −e 2 /(12πm * eff ) which characterizes the orbital response of a free electron gas. Indeed, this seems reasonable in the limit of α R ∆ xc with the chemical potential positioned in the majority band, as the true orbital magnetic susceptibility of the Rashba model (as calculated by Fukuyama's formula 20,21 ) reduces to χ oms ∼ χ ↑+↓ LP /2 in the same limit. For |µ| ≈ |∆ xc | we therefore arrive at at the intuitive result guided by Eq. (2) with B eff replaced by the chiral emergent field: This reflects the fact that in the limit of |α R | |∆ xc | the emergence of chiral orbital magnetization can be understood as the coupling of a mixed SU (2) gauge field to the diamagnetic Landau-Peierls susceptibility.
The behavior of COM becomes complicated and deviates remarkably from the α R -linear expression given by Eq. (12) as the Rashba parameter increases. To demonstrate this, we numerically calculate the value of M com at the center of a skyrmion, in a wide range of parameters ∆ xc and α R of the Rashba Hamiltonian, Eq. (5), while fixing the chemical potential at µ = 0. We parameterize the skyrmion in the polar coordinates (ρ, φ) by choosinĝ n(ρ, φ) = (sin θ(ρ) cos Φ(φ), sin θ(ρ) sin Φ(φ), cos θ(ρ)) T as the local magnetization vector 1 . Here, we define Φ(φ) = mφ + γ with the vorticity m and the helicity γ. For a Néel skyrmion γ = 0, whereas a Bloch skyrmion is represented by the value γ = π/2. The topological charge of the skyrmion then equals N sk = dxdyn · (∂ xn × ∂ yn )/(4π) = −m. In order to model the radial dependency refer to Romming et al. 33 and choose a 360 • domain wall, which is described by two parameters: the domain wall width w and the core size c (see Methods).
The results are presented in Fig. 2 for a Néel skyrmion (γ = 0) with w = 20 nm, c = 0 nm and m = 1. The magnetization is given in units of µ * B /nm 2 with the effective Bohr magneton µ * B = e /(2m * eff ). In this plot, we observe that while the gauge field picture is valid in the limit of ∆ xc /α R → ∞, there exists a pronounced region in the (α R ,∆ xc )-phase-space where COM exhibits a strong non-linear enhancement. This is in contrast to the case of Bloch skyrmions, where COM vanishes identically for all (α R , ∆ xc ), reflecting the symmetry of the Rashba coupling. It also elucidates our terminology, since already the gauge field description can be used to verify that M com ∝ cos γ (for vorticity m = 1), thereby making COM explicitly dependent on the helicity.
Topological Orbital Magnetization. The TOM appears as the correction to the OM which is second order in the gradients of the texture, Eq. (4), and while it vanishes for onedimensional spin-textures, we show that it is finite for 2D textures such as magnetic skyrmions. In contrast to COM, the TOM is non-vanishing even without spin-orbit interaction. To investigate this, we set α R to zero, reducing the effective vector potential to A = A xc and with the emergent field turning into (B xc eff ) z , Eq. (9). The gradient expansion (see Methods) Figure 4. a) When αR ∆xc, the integrated mtom is a topological quantity which, in particular, cannot distinguish between topologically equivalent structures such as Néel (black triangles) and Bloch skyrmions (black stars). Tracing the vorticities m = 1, 2, 3, 4 as a function of αR and with ∆xc = 0.9, µ = 0, the top figure illustrates how mtom (in units of m0 = −µ * B /12) passes from its regime of topological quantization (mtom/m0 = −m) to a regime of strong enhancement, with Néel and Bloch structures clearly separated. Intermediate phases form a continuum between the two values (shaded regions). The inset demonstrates for the case of Néel skyrmions, that a level structure is still present at αR = 2 eV. b) For the particular case of αR = 2.0 and ∆xc = 0.9, µ = 0 the γ phase shift is used to interpolate from a Néel to a Bloch type Skyrmion, leading to a drastic loss of mtom. Variations in the shape of the skyrmion, as quantified by the ratio of c/w, have a very small effect. now indeed reveals that which again confirms the gauge-theoretical expectation. Remarkably, the similarity between Eqs. (12) and (13) underlines the common origin of the COM and TOM in the "effective" magnetic field in the system, generated by a combination of a gradient ofn along x with spin-orbit interaction (in case of COM), and by a combination of the gradients ofn along x and y (in case of TOM).
To explore the behavior of TOM in the presence of spinorbit interaction, α R = 0, we numerically compute the value of TOM at the center of the Néel (Bloch) skyrmion with para-meters used in the previous section, as function of ∆ xc and α R (at µ = 0). The corresponding phase diagram, presented in Fig. 2, displays two notable features. The first one is the relative stability of Eq. (13) against a perturbation by a spin-orbit field in the limit of |∆ xc | |α R |. The second one is the significant enhancement of TOM in the regime where |α R | > |∆ xc |, similar to COM (albeit over a larger part of the parameter space).
Interplay of topologies. The phase diagrams in Fig. 2 have been evaluated at the core of a skyrmion. We now take a more global perspective and analyze the decomposition of the overall orbital magnetization into its constituent parts M com and M tom as a function of ρ, the radial position inside the skyrmion with w = 20 nm and c = 0 (see Fig. 3b)). By fixing α R to 2 eVÅ, and ∆ xc to 0.9 eV with µ = 0, we position ourselves precisely in the region of orbital enhancement discussed above in the context of the phase diagrams. When the local direction of the magnetization (parametrized by spherical coordinates θ and φ) is close to the z-axis, M com and M tom are rather gently varying whereas their behaviour reveals strong resonance in the vicinity of in-plane directions (θ ≈ π/2). The emergence of this resonance coincides with an occurrence of a band crossing at the critical k-value of k c = |∆ xc /α R | with the polar coordinate in the Brillouin zone of φ k = φ − π/2 in the local ferromagnetic electronic structure which corresponds to the given magnetization direction, see Fig. 3a).
It is known that this specific band crossing in the Rashba model leads to a vastly enhanced diamagnetic susceptibility 27 and in close analogy, a strong response in M com and M tom can be expected based on Eq. 2. To study the origin of this effect in greater detail, we plot M tom as a function of µ for two different magnetization directions. The results, presented in Fig. 3d) and e) reveal the sensitivity of M tom to the SOImediated deformation of the purely parabolic free-electron bands separated by ∆ xc . The magnitude of TOM is largest and exhibits pronounced oscillations in a narrow energy interval around the band edges. When we turnn into the in-plane direction, it can be seen how the resonances of M tom are enhanced in magnitude and are carried along by those band extrema which eventually touch at θ = π/2, pushing the peaks of M tom through the chemical potential (which was aligned to µ = 0 for Fig. 3b)). For three different values of the chemical potential (indicated by the symbolic arrows) the strongly µ-dependent real-space density of the total orbital magnetization M is shown in Fig. 3f). This anisotropic behaviour cannot be accounted for within the emergent magnetic field picture which only relies on the real-space texture with its associated topological charge and winding density.
The "critical" metallic point in the Rashba model that we encounter is topologically non-trivial. Indeed, the upper and lower bands of the magnetic Rashba model bare non-zero Chern numbers, C 1 = ±sgn(∆ xc )1/2, with the sign depending on the band 34 . The Chern number is a topological invariant of energy bands in k-space and can only change when bands are crossing. Since the sign of C 1 changes under the transformation ∆ xc → −∆ xc , the emergence of the critical metallic point at θ = π/2 is enforced when the direction of the magnetization is changed from θ = 0 to θ = π. This is illustrated in Fig. 3c). In the context of topological metals, such a point is known as a mixed Weyl point 35 , owing to the quantized flux of the Berry curvature permeating through the mixed space of k and θ (as confirmed explicitly by the calculations for the magnetic Rashba model). These points have recently been shown to give rise to an enhancement of spinorbit torques and Dzyaloshinskii-Moriya interaction in ferromagnets 35 . Here, we demonstrate the crucial role that such topological features in the electronic structure could play for the pronounced chirality-driven orbital magnetism of spin textures. Given the observation that TOM simply follows the evolution of the electronic structure in real space via the direction of the local magnetization (as illustrated in a schematic way in Fig. 3d)), the close correlation of real-and reciprocal space topologies offer promising design opportunities in skyrmions or domain walls of transition-metals with complex anisotropic electronic structure.
Topological Quantization and Stability. One of the key properties of M tom is its origin in the local real-space geometry of the texture. This has drastic consequences for the topological properties of the overall orbital moment of twodimensional topologically non-trivial spin textures as we reveal below for the case of chiral magnetic skyrmions. We thus turn to the discussion of the total integrated values of the orbital moments in chiral spin textures by defining them as Concerning the total value of the COM-driven orbital moment in one-dimensional uniform 360 • or 180 • chiral domain walls it always vanishes identically by arguments of symmetry (although it can be finite for example in a 90 • wall). In sharp contrast, the TOM-driven total orbital moment of isolated skyrmions generally does not vanish. This can be most clearly shown in the limit when the gauge-field approach is valid (i.e., ∆ xc α R ). In this case, the total flux of the emergent topological and chiral fields through an isolated skyrmion is given by It then follows from Eq. (12), that the integrated value of M com vanishes while Eq. (13) predicts the quantization of the topological orbital moment m tom to integer multiples of χ ↑+↓ LP Φ 0 = −µ * B /6 (at |µ| = |∆ xc |). In this limit, the skyrmion of non-zero topological charge N sk = 0 thus behaves as an ensemble of N sk effective particles which occupy a macroscopic atomic orbital with an associated value of the orbital angular momentum of µ * B /6. This quantization is explicitly confirmed in Fig. 4a), where we present the calculations of m tom for Néel and Bloch-type skyrmions with dimensions c = 0 nm and w = 20 nm at a fixed value of ∆ xc = 0.9 eV while varying α R for different topological charges N sk ∈ {−1, −2, −3, −4}. The results, presented in units of m 0 = −µ * B /12 (corresponding to the Landau-Peierls limit at µ = 0 and N sk = −1), reveal a stable plateau, corresponding to the regime of topological quantization, where m tom attains the value N sk m 0 .
In the opposite limit of α R > ∆ xc the magnitude of M tom can be enhanced drastically with respect to the topologically quantized value. When α R reaches a magnitude of about 1 eVÅ, the emergent field picture breaks down and we discover a drastic increase in Néel-m tom by as much as one order of magnitude upon increasing α R . And although m tom is not topologically quantized in this regime, it still attains a distinctly different magnitude for different skyrmion charges, and while it is weakly dependent on the c/w ratio (up to a couple of percent), when c = const the variations of w keep m tom strictly constant (see the insets in Fig. 4 with c = 0). The latter robustness can be demonstrated already analytically on the level of Eq. (4) using the transformation of coordinates x → x/w. Remarkably, in the regime of enhanced SOI, the strong dependence of the local TOM on the helicity of the skyrmion (i.e. Néel or Bloch), uncovered in Fig. 2, is translated into a drastic dependence of the overall topological orbital moment on the way that the magnetization rotates from the core towards the outside region, as shown in Fig. 4. Such behavior of the topological orbital moment with respect to deformations of the underlying texture suggests that monitoring the dynamics of the orbital magnetization in skyrmionic systems can be used not only to detect the formation of skyrmions with different charge, but also to distinguish various types of dynamical "breathing" modes of skyrmion dynamics 36 .

Discussion
On a fundamental level, COM and TOM arise as a consequence of the changes in the local electronic structure caused by a non-collinear magnetization texture. Since the effective magnetic fields directly couple to the orbital degree of freedom, they lead to the emergence of chiral and topological orbital magnetization. While this intuitive interpretation in terms of real-space gauge fields eventually breaks down at large SOI, it makes room for a regime of strong enhancement in which the intertwined topologies of real-and reciprocal space lead to novel design aspects in the bandstructure engineering of orbital physics. This is possible by exploiting either the spin-orbit interaction or the dispersive behaviour of the bands, i.e. their effective mass. In particular, the metallic point in the mixed parameter space of the non-collinear Rashba model reveals its strong impact on COM and TOM. Such critical points will have a pronounced effect on the orbital magnetism even if they are emerging on the background of metallic bands in transition-metal systems. Our analysis therefore indicates in which materials an experimental detec-tion of orbital magnetism associated with the non-collinearity is the most feasible. By numerically evaluating the magnitude and real-space behavior of the TOM and COM, we thereby uncover that by tuning the parameters of surface and interfacial systems the orbital magnetism of domain walls and chiral skyrmions can be engineered in a desired way.
Concerning experimental observation of the effects discussed here, M tom and M com could be accessible by techniques such as off-axis electron holography 37 (sensitive to local distribution of magnetic moments), or scanning tunneling spectroscopy (sensitive to the local electronic structure) in terms of B-field induced changes in the dI/dU or d 2 I/dU 2 spectra 38 . An already existing proposal for the detection of non-collinearity driven orbital magnetization of skyrmions by Dias et al. relies on X-ray magnetic circular dichroism (XMCD) which is able to distinguish orbital contributions to the magnetization from the spin contributions 9 .
Further, the emergence of COM and TOM can give a thrust to the field of electron vortex beam microscopy 39 -where a beam of incident electrons intrinsically carries orbital angular momentum interacting with the magnetic system -into the realm of chiral magnetic systems. For example, we speculate that at sufficient intensities, electron vortex beams could imprint skyrmionic textures possibly by partially transforming its orbital angular momentum into TOM. Since the topological orbital moment is directly proportional to the topological charge of the skyrmions, we also suggest that the interaction of TOM with external magnetic fields could be used to trigger the formation of skyrmions with large topological charge. Ultimately, the currents of skyrmions can be employed for lowdissipation transport of the associated topological orbital momenta over large distances in skyrmionic devices.
While in this work we focus primarily on TOM, the here discovered chiral orbital magnetization has been an overlooked quantity in chiral magnetism so far. Besides the fact that it emerges already in one-dimensional chiral systems and serves as a playground to study the effects of mixed space Berry phases, it can reach very large values depending on the details of the texture as well as strength of SOI. Even in case of skyrmions, where the argument of vanishing effective flux, Eq. 16, might suggests that COM is not of importance, it turns out that beyond the α R → 0 limit the integral effect of M com can be substantially enhanced in a way similar to TOM. A prominent example for the importance of COM is given e.g. in Vanadium-doped BiTeI 40 which has a large SOC of α R = 3.8 eVÅ and an exchange gap of ∆ xc = 45 meV. If this material would host Néel skyrmions (m = 1, w = 20 nm, c = 0 nm), m com would reach approximately 12µ * B which is magnitudes larger than the corresponding m tom of about −0.7µ * B . Creating skyrmions and large COM in strong Rashba systems might therefore be a promising direction to pursue.
In a wider perspective, the emergence of TOM and COM gives rise to a physical object which is directly connected to the orbital degree of freedom with the advantage that it can be understood from a semiclassical perspective in a way which is engineerable and controllable. Our findings thus open new vistas for exploiting the orbital magnetism in chiral magnetic systems, thereby opening interesting prospects for the field of "chiral" spintronics and orbitronics.

Methods
Gradient expansion. The expansion in exchange field gradients is naturally achieved within the phase-space formulation of quantum mechanics, the Wigner representation 17,22 . The key quantity in this approach is the retarded single-particle Green's function G R , implicitly given by the Hamiltonian H via the Dyson equation where x µ = (t, x) and π µ = ( , π) are the four-vectors of position and canonical momentum, respectively. The latter of the two, in terms of the elementary charge e > 0 and the electromagnetic vector potential A, is related to the zero-field momentum p by the relation πµ(x, p) = pµ + eAµ(x). The -product, formally defined by the operator of left-and right-acting derivatives ↔ ∂ , allows for an expansion of G R in powers of , gradients ofn and external electromagnetic fields, captured in a covariant way by the field tensor F µν = ∂x µ A ν − ∂x ν A µ22 .
In this work, we are after the orbital magnetization (OM) in zdirection. Given the grand canonical potential Ω, the surface density of the orbital moment is given by 25,41 M (x) = −∂B Ω(x) , which requires an expansion of Ω up to at least first order in the magnetic field B = Bez in the collinear case. In the limit of T → 0, the grand potential is asymptotically related to the Green's function G R via where denotes the imaginary part, the integral measure is defined as dp = d d 2 p, f ( ) represents the Fermi function f ( ) = (e β( −µ) + 1) −1 , µ is the chemical potential and β −1 = kBT . In our approach, deviations from the collinear theory enter the formalism as gradients ofn and can be traced systematically in G R and in Ω , finally leading to Eq. (3) and (4).
Computational details. All calculations were performed with a Green's function broadening i0 + → iΓ with Γ = 100 meV while we approach the zero-temperature limit by setting kBT = 10 meV. The k-space integrals are then performed on a quadratic 512 × 512 mesh. The effective electron mass was set to m * eff /me = 3.81 everywhere except for the example of V-doped BiTeI, where m * eff /me = 0.1 42 .
The topological charge is then given by N sk = 1 4π dx dyn · (∂xn × ∂yn) Assuming Φ(φ) = mφ+γ, with the vorticity m ∈ Z and the helicity γ ∈ R (Néel skyrmions correspond to γ = 0, Bloch skyrmions to γ = π/2), as well as the property θ(0) = π and θ(∞) = 0, the integral evaluates to N sk = −m. A realistic profile satisfying these requirements and which is used in this work is given by 33 with the core size c and the domain wall width w.
Data availability. The code and the data that support the findings of this work are available from the corresponding authors on request.