Mewen Crespo
email: mewen.crespo@univ-rennes.fr

Guy Casale
email: guy.casale@univ-rennes.fr

Loïc Le Marrec
email: loic.lemarrec@univ-rennes.fr

Continuum mechanics of defective media: an approach using fiber bundles ⋆

Keywords: Generalized continuum, Fiber bundles, Connection, mechanics

The kinematics of a micro-structured material is geometrically modeled through the framework of fiber bundle geometry. The material continuum is a fiber bundle M B where B is compact and orientable. It is commonly agreed that connections with curvature and torsion can describe defect densities in micro-structured materials. The aim of this work is to introduce a method to derive these objects from the kinematics in an intrinsic way. The material bundle M is therefore placed in the Euclidean fiber bundle E ≡ TE E using a placement map φ : M E. A first-order transformation F : TM TE generalizing Tφ is then introduced. Finally, using F, a metric on B, a connection on M and a solder form on M are inferred from the Euclidean structure on E. These new objects are grouped into a single one, called a pseudo-metric, which allows us to describe the current state of matter through, among other things, the curvature (disclinations) and torsion (dislocations) tensors. On one hand, we see that the torsion tensor can be non-zero even in the holonomic F = Tφ case. On the other-hand, in order for the material to have a non-zero curvature tensor, we see that one must have a non-holonomic first-order transformation: that is, F ̸ = Tφ.

Introduction

We are interested in the study of micro-structured materials. That is, materials in which the microscopic state at a point plays an import role in its macroscopic behavior and therefore have to be included in its geometric description. Mathematically, an n-dimensional micro-structured material with a k-dimensional micro-structure (or a n × k material for short) is therefore a fiber bundle manifold M π M -B of rank k over a compact orientable n-dimensional manifold B, ⋆ Supported by the ANR-11-LABX-0020-0 program Henri Lebesgue Center.

which is the classical macroscopic body.

At a point b in the macroscopic body B, the micro-structure of a microelement is represented by the fiber M b of M over b. Those elements are interpreted as first order infinitesimal neighborhoods of geometrical macroscopic points (eventually embedded in a bigger space of higher dimension), as it is commonly considered for crystal modeling. For this reason, M is often closely related to a tangent space (not necessarily TB). In our case, the material is embedded in the Euclidean space. Therefore, n ∈ {1, 2, 3} and k = 3.

Fig. 1: The standard horizontal representation of a generic material bundle. The macroscopic space is represented horizontally and micro-elements are shown to be inside the macroscopic points. Note that neither B nor M b is assumed to be 1-dimensional, this is only done for illustration purpose.

Examples of micro-structured materials include the 1-dimensional beams (1×3 materials), the 2-dimensional plates (2×3 materials) and the 3-dimensional Cosserat continuum (a 3 × 3 material). Some authors [5, p138-141] choose to use the principal bundle (see section 2.2). It is a fiber bundle over the same base manifold B whose fibers are not the micro-elements but their possible configurations (a sub-group of GL( 3)). This is why, at first, those examples may seem to have different micro-structure as, in practice, one often ends up manipulating the principal bundle. Those could for example, depending on assumptions made on the material, be isomorphic to the 1-dimensional circle SO(2) ≃ S 1 (for planar inextensible beams), the 2-dimensional space SO(2) × R * (for planar extensible beams) or the 3-dimensional space SO(3) (for some plates and Cosserat continuum).

The standard interpretation

Figure 1 shows the standard representation of a generic material bundle. This interpretation of the mathematical model is central to our problem. It states that micro-elements are inside macroscopic points and are therefore imperceptible to an experimenter. In general, this interpretation states that the projection, to the macroscopic spaces B and E, of any kinematic object of our model can be seen as a kinematic object of the classical continuum mechanics. For example, this implies that φ : B E (defined later on) corresponds to an usual macroscopic placement map that can be observed and measured by an experimenter.

That the macroscopic kinematics is classical does not mean that the macroscopic dynamics will be. Macroscopic kinematic object are independent of the micro-structure, but microscopic or mixed objects can depend on the macroscopic states. This means that a coupling between the macroscopic and microscopic variables can append in the dynamics' equations. In other words, in general, φ will not behave classically.

1 First formalism

Placement map

By analogy with M, the ambient space is modeled as the Euclidean fiber bundle E ≡ TE π E -E of rank 3 over the usual ambient space E, the 3-dimensional Euclidean space. We chose a Lagrangian approach. Accordingly, one defines a placement map. This latter is a fiber bundle embedding φ : M E which maps the micro-element over b ∈ B to the micro-place over φ(b) ∈ E where φ : B E is called the shadow of φ (see the commutative diagram bellow). This placement map φ defines a macroscopic material placement map, via its shadow φ and the differential Tφ.

TE TE E E TE TM TE B E TM Tπ E π TE π E π TE Tπ M F π TB Tφ φ π M φ π TM F
In order to also have a description of the material's defects, one must have a generalized material placement map. To complete the description of a generalized material placement map, one additionally has to introduce a generalized first order placement map F : TM TE. The standard interpretation then requires F (VM) ⊂ VE. Physically, this means that a microscopic change in the material causes a microscopic change in the ambient image. In particular, this implies that F can be seen as a fiber bundle morphism between TM The standard interpretation then also requires that the shadows of F with respect to the usual structures and these new structures are, respectively, φ : M E and Tφ : TB TE. Those two requirements reduces the degrees of freedom of F to k • (k + n) ∈ {12, 15, 18} real coefficients, which are scalar fields on M (i.e. F lives in a (k +n)×(k •(k +n)) vectorial fiber bundle). In the classical case F ought to be the differential of φ -i.e. F = Tφ -but as it turns out, this assumption is too restrictive and corresponds to the absence of disclinations in the micro-structure (see [START_REF] Nguyen | On tangent geometry and generalized continuum with defects[END_REF]).

Metric, connection and solder form

The existence of a first order placement map allows one to superimpose the tangent of the body over the tangent of the ambient space. Doing so, first order functions of the ambient space such as tensors, can be seen in the material space. Mathematically, this correspond to a pull-back by the first order placement map.

First, there exists a metric g on E: the Euclidean metric. This allows the measure of angles and lengths. Using the first order macroscopic placement Tφ, one can pull it back as a metric G on M, which gives the Cauchy-Green tensor (see 1 the equation and diagram below). In the case of generalized continua, the macroscopic metric is accompanied by two micro-structured linear applications: the connection and the solder form.

∀b ∈ B G b := T * φ(b) φ • g φ(b) • T b φ T * B T * E TB TE T * φ G Tφ g
Figure 2 shows the vertical representation of a generic connection on a material bundle. A connection on E can be defined as a right-inverse Γ : TB TM of Tπ M : TM TB (i.e. Tπ M • Γ = Id TB -please note that other equivalent definitions also exist [START_REF] Epstein | The Geometrical Language of Continuum Mechanics[END_REF][START_REF] Michor | Topics in Differential Geometry[END_REF][START_REF] Marsden | Mathematical Foundations of Elasticity[END_REF]). While it does not induce a notion of length or angle, it still induces a notion of local geodesics and, as shown on the figure, parallel transport of micro-elements along the horizontal lift of a macroscopic path. This later, allows to connect close microscopic spaces by following the flow generated by Γ. If the macroscopic path is a closed curve in B, one may witness some lack of closure of the horizontal lift in M. This lack of closure is measured by the curvature and torsion in a way similar to the way Frank and Burgers vectors are introduced and used in crystallography (see [START_REF] Katanaev | Geometric theory of defects[END_REF]).

On E, one has the Levi-Civita connexion, which is defined as the unique connexion with no torsion whose geodesics are those of the Euclidean metric. Coordinate-wise, it corresponds to the inclusion TE ≃ E E × {0} ⊂ E × R 6 ≃ TE. In a similar way to what has been done for the metric, the connection on M, 1 Here T * φ : T * E T * B is the transpose of the differential Tφ : TB TE of φ : B E.

which is not a tensor, is then obtained by pulling back this Euclidean Levi-Civita connection γ of E into a connection Γ on M. This involves the use of the fiber bundle morphism F : TM TE to superimpose TM over TE. The formula and associated diagram are as follows:

∀b ∈ B, ∀m ∈ M b , Γ m := F b -1 • γ φ(m) • T b φ TM TE TB TE F Tπ M Tπ E Γ Tφ γ
A solder form on E is, in an analog fashion, an injective2 map ϑ : T x E

V p E for any x ∈ E, p ∈ E x . While a connection glues microscopic grains together, a solder form glues a grain to the macroscopic space. From a vector at a macroscopic point, it generates a vector field inside the microscopic grain at that macroscopic point. As an ambient solder form, we choose the canonical Euclidean solder form ϑ which comes from the identification E ≃ TE. Just like for γ, ϑ can be pulled back onto M as Θ : TB VM. The formula and associated diagram are as follows:

∀b ∈ B, ∀m ∈ M b , Θ m := F b -1 • ϑ φ(m) • T b φ VM VE TB TE F ϑ -1 Θ Tφ ϑ
Fig. 2: The vertical representation of a generic material bundle. The macroscopic space is represented horizontally and micro-elements are unfolded vertically above it. Note that neither B nor M b is assumed to be 1-dimensional, this is only done for illustration purpose.

Pseudo-metric

Another way to interpret the connection and solder form is as embeddings. The connection can be seen as an application which takes a macroscopic vector v ∈ TE and embeds it in the total tangent space as γ (v) ∈ TE which can be interpreted as the corresponding vector with the same macroscopic part

(Tπ E • γ • v = v) and no vertical 3 part (v γ • γ • v = 0).
Similarly, ϑ(v) ∈ VE ⊂ TE can be interpreted as the corresponding vector with the same vertical part ϑ

-1 • v γ • ϑ • v = v and no horizontal part (h γ • ϑ • v = 0). In this interpretation, ϑ -1 • v γ : TE
TE can be seen as a "microscopic" projection. In particular,

Tϖ = Tπ + ϑ -1 • v γ : TE TE
is also a projection which will be called the projection of interpretation as it gives the sum of the macroscopic part and the microscopic part, both seen in the macroscopic space TB, which is how one interprets v.

This projection of interpretation gives an alternative approach for the generalization of the metric. Indeed, one can simply pull-back the Euclidean metric g : TE T * E onto TE using Tϖ. One then obtains

g : TE -T * E := T * ϖ • g • Tϖ
Its kernel is the kernel of Tϖ. Those are the vectors whose macroscopic part is the opposite of their microscopic part. In particular their microscopic and macroscopic part are of the same magnitude, which mean these are not physically acceptable vectors. The kernel of g, as it turns out, is the horizontal space H γ-ϑ = Im (γ -ϑ) of the connection γ -ϑ. The solder form ϑ is, by construction, isometric with respect to g | VE : VE

V ⋆ E and g. This means that, from g, one can extract γ -ϑ and ϑ up to a "rotation" in O g, g | VE .

The tensor g is called the ambient pseudo-metric as it is symmetric and positive but its kernel is non-trivial (otherwise it would be a metric). As g with Tφ, g can be pulled-back on TM using F. One obtains the material pseudometric

G : TM -T * M := F * • g • F
Similarly to the ambient case, one has that ker (G) = H Γ-Θ and that Θ is an isometry for G and G | VE . This means that, if G is known, one can extract the connection Γ -Θ and the orbit of

Θ under O G, G | VE from G. Nevertheless, G : TB T * B cannot be extracted from G : TM T * M directly, unless
F : TM TE is entirely known (in which case G is also known). The pseudometric G is a canonical extension of the right Cauchy-Green tensor G. As G plays a central role in the classical case, one expect G to also play a key role in the generalized case. These observations therefore justify the use of G, Γ and Θ as our primary tools.

2 Constitutive equations and structure group

Stress and constitutive equations

By formulating continuum mechanics using differential geometry (as in [START_REF] Kolev | Éléments de géométrie différentielle à l'usage des mécaniciens[END_REF][START_REF] Epstein | The Geometrical Language of Continuum Mechanics[END_REF][START_REF] Gonzalez | A First Course in Countinuum Mechanics[END_REF]) one obtains a geometric formulation of static equilibrium that can be generalized to the micro-structured case. For example, the macroscopic Cauchy stress S and the macroscopic second Piola tensor Σ are:

S : T * E -TE Σ : T * B -TB
In the generalized version, tensors on E (resp. B) naturally become tensors on E (resp. M). This gives

S : T * E -TE Σ : T * M -TM
Which, for consistency, are required to be morphisms for both structures TE -TB . This is a manifestation of the standard interpretation, just like it was for F, which also requires that the shadows be Id E and S resp. Id M and Σ .

The notion of elasticity and hyper-elasticity can be generalized to the microstructured case. The elasticity simply states that S (resp. Σ) is a linear function S resp. Σ of F, or a linear function of G in the isotropic case. The hyperelasticity on the other side states that this function Σ is obtainable as the differential TW of a functional W of F (resp. G in the isotropic case) called the energy density. Such an energy, in the isotropic case, can depend on invariants of G, Γ and Θ. Already existing formulae for the macroscopic case can be extended by replacing G with G. For example, identifying TM and T * M for simplicity (this is roughly equivalent to having a reference pseudo-metric), the Saint Venant-Kirchhoff constitutive law gives, under the substitution G G:

W (G) := a 4 tr (G -Id) + b 4 tr (G -Id) 2 Σ := a tr (G -Id) Id + b (G -Id)
Under a specific frame induced by Γ, G splits into blocks as

Tπ * M • G • Tπ M 0 0 (Id -Γ • Tπ M ) * • G • (Id -Γ • Tπ M )
Hence, we see that even with a naive extension of the Saint Venant-Kirchhoff formula, there is a micro/macro coupling occurring through the term

tr (G -Id) = tr (G) -k + tr (Id -Γ • Tπ M ) * • G • (Id -Γ • Tπ M ) -3
If we place ourselves in the holonomic case φ ≡ Tφ and F = Tφ. That is, if the microscopic kinematics is a copy of the macroscopic one. Further computations show that one has tr(G) = 2tr(G) and, up to a constant rescaling of a and b, Σ is just two (diagonal) copies of Σ. That is, the microscopic part of the stress is just an independent copy of its macroscopic part and therefore the dynamics is classical. In general, however, it is not the case.

It should be noted that, simply extending classical formulae this way does not provide all possible energies, as other invariants exist. For example, even tr (G) will not be accessible as it will be replaced with tr(G) in the generalized case. Furthermore, this would not necessarily lead to geometrically exact laws as the micro-structure is ignored. Therefore, invariants of the pseudo-metric G need to be used instead. That is, invariants of the connection Γ and the torsion Θ on M must be used, in addition to those of the metric G on B. In theory, this could allow the formulation of geometrically exact constitutive equations making use of the micro-structure. Equations which could use, for instance, the torsion and curvature tensors of the connection to quantify these generalized strain and stress.

Structure group and principal formalism

Although we took care to omit this detail until now, fiber bundles are not merely spaces locally isomorphic to Cartesian products. They are actually equipped with a group, called the structure group [START_REF] Epstein | The Geometrical Language of Continuum Mechanics[END_REF], acting on its typical fiber. This group can be physically interpreted as the group of admissible change of coordinates or the set of microscopic configurations. Objects of the theory are therefore required to interact properly with it. The structure group is usually simpler than the typical fiber of M. For this reason, using the principal bundle associated to M is sometimes easier. The principal bundle associated to a fiber bundle M B is a fiber bundle P M whose base space, structure group and (in some sens) transition maps are the same but whose typical fiber is the structure group G M itself. For example, in the case of planar and inextensible Timoshenko-Ehrenfest beams, the structure group is SO(2) which is isomorphic to the circle S 1 . Therefore, the principal material bundle is P M B whose typical fiber is G M ≃ S 1 .

Objects defined on M and E can canonically be defined (or transferred) on P M and P E (see [START_REF] Epstein | The Geometrical Language of Continuum Mechanics[END_REF][START_REF] Marsden | Mathematical Foundations of Elasticity[END_REF]). A placement map φ : M E can for example be decomposed as a macroscopic placement map φ : B E and a material section (i.e. a microscopic configuration) µ : B G M . In the case of the planar and inextensible Timoshenko-Ehrenfest beam, we see that a material placement map is therefore described by a classical placement map φ : B E and an angle θ b ∈ S 1 ≃ [0, 2π[ at each b ∈ B (see figure 3).

-

  TE (notice the change of base spaces and projection maps, see the commutative diagram above).

-

  TE resp. TM π TM -M and TM Tπ M

Fig. 3 :

 3 Fig. 3: Deformation of a Timoshenko beam. The normal section rotates by an amount θ x with ∂w ∂x ̸ = sin (θ x ) ∂w ∂x is the z-aligned part of Tφ . Source: Wikipedia

A quick comparison of dimensions gives dim(TxE) = 2 dim(E) = 6 = 2 rank(E) = dim(VpE). Hence, in the Euclidean case, ϑ is a bijection.

We use here the standard notation hγ = γ • TπE and vγ = IdTE -hγ such that TπE • hγ = TπE and TπE • vγ = 0.

Acknowledgements We would like to thank the C.N.R.S' research group GDR GDM C.N.R.S. for stimulating the interactions between differential geometry and mechanics. An interaction which is at the heart of this work. We would also like to thank CNRS, Univ Rennes, ANR-11-LABX-0020-0 program Henri Lebesgue Center for their support.


Orbital Hall physics in two-dimensional Dirac materials Armando Pezo, * Diego García Ovalle, and Aurélien Manchon † Aix-Marseille Université, CNRS, CINaM, Marseille, France.

(Dated: December 21, 2023)

Orbitronics has recently emerged as a very active research topic after several proposals aiming to exploit the orbital degree of freedom for charge-free electronics. In this communication, we investigate orbital transport in selected two-dimensional systems to better understand which parameters govern the intra-atomic and inter-atomic contributions to the orbital Hall effect. We study the impact of the gap, the role of the materials' topology and the influence of the disorder on spin and orbital Hall transport. Starting from the Kane-Mele model, we describe how the orbital moment behaves depending on the material's topology and clarify the influence of the gap on the orbital Hall conductivity. We then extend the study to realistic topologically trivial and non-trivial materials, and find that the topology has little qualitative influence on the orbital Hall conductivity. In contrast, we observe that the energy dispersion has a more dramatic impact, especially in the presence of disorder. Remarkably, our results suggest that the intra-atomic orbital Hall current is more robust against scattering than the inter-atomic one, without further impact of the topological properties of the system under consideration.

I. INTRODUCTION

Recent theoretical and experimental efforts suggest that the orbital angular momentum of electrons can be used as an alternative degree of freedom to the spin angular momentum [1][2][3][4][5][6]. In contrast with the generation of spin currents that necessitates either a ferromagnet or a heavy metal, orbital currents can be induced electrically using light metals, thereby presenting a potential technical advantage in terms of materials scarcity [1,3]. Current research is being developed along two directions. A first direction takes advantage of the vast experience acquired on spin transport in transition metal heterostructures [7,8]. Tight-binding and first principles calculations have suggested that certain light metals such as V, Cr or Cu can host large orbital Hall effect [5,9], resulting in the experimental demonstration of orbital torque and magnetoresistance [10][11][12][13]. In the last years new experimental developments have unlocked the synthesis of two-dimensional materials opening new possibilities for spintronics [14]. To date, most of the attention has been focused on graphene and transition metal dichalcogenides, where valley Hall effect and orbital Hall effect coexist [15][16][START_REF] Cysne | Rappoport[END_REF][18].

In most early theoretical studies on the orbital Hall effect, the orbital moment was assumed to be mostly of intra-atomic origin, adopting the so-called atom-centered approximation (ACA) [1,2,5,9,[START_REF] Cysne | Rappoport[END_REF]19]. However, important developments in the theory of orbital magnetism have demonstrated that ACA is not sufficient to properly describe the orbital motion of quasiparticles in solids and that the inter-atomic contribution cannot be neglected [20][21][22][23]. Whereas the inter-atomic orbital moment is small in the case of bulk transition metals [23], it is particularly significant in materials like graphene where the intra-atomic orbital character of the conduction electrons vanishes [24]. Extending the "modern theory" of orbital magnetization to the orbital Hall effect, it recently appeared that the inter-atomic contribution cannot be neglected in general [6,16,18].

Building on our previous work on the modern theory of orbital Hall effect in realistic materials [6], we investigate orbital transport in selected two-dimensional systems to better understand which parameters govern the intra-atomic and inter-atomic contributions to the orbital Hall effect in these systems. In particular, we investigate the impact of the gap, the role of the materials topology and the influence of the disorder on spin and orbital Hall transport. We find that whereas the topology has little influence on the orbital Hall effect itself, the orbital transport exhibits markedly distinct behavior depending on the nature of the gap (spin-orbit coupling, staggered potential) and systematically increases upon reducing the gap size. Our results also show that the intra-atomic and inter-atomic orbital Hall conductivities experience different robustness against disorder. Whereas the inter-atomic orbital Hall contribution systematically decreases upon increasing impurity scattering, the intra-atomic contribution remains mostly unaffected in the gap region and decreases in the metallic regime, suggesting that overall the intra-atomic contribution is more robust against disorder than the inter-atomic one.

This work is organized as follows. In Section II, we briefly remind the concepts of the modern theory of the orbital Hall effect and use the two-dimensional Kane-Mele model to determine the influence of the lattice topology on the spin and orbital Hall transport. In Section III, we investigate the spin and orbital Hall effects in selected two-dimensional lattices computed using density functional theory (DFT) and show that orbital Hall transport is substantially influenced by the proximity to the gap and by the nature of the energy dispersion. Then, in Section IV, we investigate the impact of the Anderson-type disorder on the orbital Hall effect and show that intra-atomic and inter-atomic contributions behave differently. Conclusions and perspectives are given in Section V.

II. THEORY AND CONCEPTS

A. Modern theory of the orbital Hall effect

Let us remind the concepts related to the modern theory of orbital magnetization and its extension to the orbital Hall effect. In crystals, the orbital motion of an electron arises from the self-rotation of the wavepacket in the unit cell which includes both the intra-atomic and interatomic contributions as mentioned above. In equilibrium, this self-rotation gives rise to the orbital magnetization as long as the time-reversal symmetry is broken [20,22]. In contrast, for materials with time-reversal symmetry, an orbital magnetization can be generated out of equilibrium as long as inversion symmetry is broken [25]. When both time-reversal symmetry and inversion symmetry are preserved though, no orbital magnetization can be induced, and only the orbital Hall effect survives.

Let us start with the real space definition of the orbital moment, L = (r × pp ×r)/4. Under the parallel transport gauge condition applied on non-degenerate bands (i.e., ⟨n| ṅ⟩ = 0), this expression can be projected on the Bloch states |u n k ⟩ and recasted in the form (see [16,26])

⟨u n k | L|u p k ⟩ = e 2g L µ B Im⟨∂ k u n k | × H k |∂ k u p k ⟩ -e 4g L µ B (ε n k + ε p k )Im⟨∂ k u n k | × |∂ k u p k ⟩. (1)
Here, |u n k ⟩ is the periodic part of the Bloch state associated with the energy ε n k , v = ℏ -1 ∂ k H k is the velocity operator, H k being the Hamiltonian in momentum space, µ B = eℏ/2m e is Bohr's magneton and g L = 1 is Landé's g-factor. This expression can be further formulated in a more tractable identity given by

⟨u n k | L|u p k ⟩ = eℏ 2 4µ B Im q̸ =n,p 1 ε q k -ε n k + 1 ε q k -ε p k ⟨u n k |v|u q k ⟩ × ⟨u q k |v|u p k ⟩ . (2)
In the Bloch state basis, the matrix element of the orbital current operator, defined as

J γ i = 1/2{L γ , v i }, reads ⟨u n k |J γ i |u m k ⟩ = 1 2 p (⟨u n k |v i |u p k ⟩⟨u p k |L γ |u m k ⟩ +⟨u n k |L γ |u p k ⟩⟨u p k |v i |u m k ⟩) , (3)
where the orbital moment L γ is either the atomic orbital moment operator (intra-atomic orbital current) or Eq. (2) (total orbital current). In the following, we will take L γ = L z since we focus on two-dimensional lattices. In the linear response theory, spin and orbital currents are time-reversal symmetric which allows one to consider the intrinsic Fermi sea contribution of the Kubo formula [27]. The orbital conductivity reads

σ z ij = -2ℏe BZ d 3 k (2π) 3 n f n (k)× Im m̸ =n ⟨u n k |J z i |u m k ⟩ ⟨u m k |v j |u n k ⟩ (ε n k -ε m k ) 2 , ( 4
)
where f n (k) is the equilibrium Fermi distribution function. The quantity that multiplies the Fermi function is usually referred to as the orbital Berry curvature, in analogy with the conventional Berry curvature where the orbital current operator is replaced by the velocity operator [3,[START_REF] Cysne | Rappoport[END_REF]28]. Due to the ubiquitous presence of spin-orbit coupling in real materials, orbital currents usually coexist with spin currents and probing orbital-only responses is experimentally challenging. To overcome this difficulty, one usually relies on materials with vanishing or low spinorbit coupling so that spin currents can be neglected [11,29]. In topological materials though, the topological transition is mostly driven by the strong spin-orbit coupling of the heavy elements, so that spin and orbital currents are entangled. A pedagogical tool to evaluate how these two effects behave in two-dimensional topological materials is the Kane-Mele Hamiltonian [30]. This model features a two-dimensional honeycomb lattice with spin-orbit coupling and its Hamiltonian reads where the creation and annihilation operators c i , c † i are spinors representing the spin degree of freedom. The first term is the usual nearest neighbour hopping (t), the second term is the spin-orbit coupling term (λ SOC ) acting on the next-nearest neighbours and the last term is the on-site staggered potential term (ϵ i ∆, ϵ i = ±1 on the different sublattices), also called Semenoff mass gap [31]. This last term leads to the appearance of an orbital magnetic moment as pointed out in Refs. [16,28]. To illustrate the connection between the Berry curvature, the orbital magnetic moment and the orbital Hall effect, let us consider the effective Hamiltonian of the Kane-Mele model, valid close to the neutrality point at K and K' points in the Brillouin zone. In the absence of Rashba interaction, the Kane-Mele model can be thought of as a double copy of the Haldane model [32], each copy corresponding to a spin sector described by the two-band Hamiltonian

H = t ⟨ij⟩ c † i c j + iλ SOC ⟨⟨ij⟩⟩ ν ij c † i s z c j + ∆ i ϵ i c † i c i , (5)
H ηs = v F (ηk x σ x + k y σ y ) + (∆ + ηs λSOC )σ z , (6)
where λSOC = 3 √ 3λ SOC , σ is the pseudospin in the sublattice space, s = ±1 refers to the spin projection and η = ±1 to the K and K' valleys, respectively. The gap is given by ∆ ηs = ∆+ηs λSOC so that the non-trivial regime is reached whenever | λSOC | > |∆|. This Hamiltonian can be written [25,33]

H ηs = Îϵ ηs k + σ • d ηs k , (7)
where ϵ ηs k is the energy dispersion of the individual bands and d ηs k describes the hybridization between bands.

Therefore the Berry curvature and orbital magnetic moment read [25] Ω ηs k,i = ±

ε ijk 2(d ηs k ) 3 d ηs k • ∂d ηs k ∂k j × ∂d ηs k ∂k k , (8)
m o,ηs k,i = - e ℏ ε ijk 2(d ηs k ) 2 d ηs k • ∂d ηs k ∂k j × ∂d ηs k ∂k k , (9)
with ± referring to the conduction and valence band, respectively and d ηs k = |d ηs k |. Equations ( 8) and ( 9) show that the Berry curvature and the orbital magnetic moment possess a very similar structure in momentum space [16,20], revealing that the orbital motion finds it origin in the Berry curvature of the Bloch states. Since the Kane-Mele model, Eq. ( 5), describes an effective four-band model taking into account the spin degree of freedom, the Bloch state does not carry any atomic orbital momentum (in other words, in this model the Bloch states are typically s or p z ). For this reason, the orbital moment comes entirely from the details of the band structure intimately correlated to the Berry curvature, as displayed as a color gradient in Fig. 1(a). Using the effective twoband Hamiltonian, one obtains

Ω ηs k,z = ± ηv 2 F ∆ ηs 2(∆ 2 ηs + v 2 F k 2 ) 3/2 , (10)
m o,ηs k,z = - e 2ℏ ηv 2 F ∆ ηs ∆ 2 ηs + v 2 F k 2 . ( 11
)
One sees that both the Berry curvature and the orbital magnetic moment are inversely proportional to the band gap [16]. The Berry curvature hot-spots appearing at different valleys characterize the topological transition: in the non-trivial phase both spin partners at the same inequivalent point (K or K') display opposite values of the orbital moment, while they have the same sign in the trivial regime, as seen in Fig. 1(b). Following the procedure outlined by Bhowal and Vignale [16], the orbital Hall conductivity for spin s and valley η reads

σ ηs OH = e 2π m e v 2 F 6ℏ 2 g L ∆ 2 ηs (∆ 2 ηs + v 2 F k 2 ) 3/2 , (12)
Therefore, the total orbital Hall conductivity is the sum of individuals contributions, σ OH = η σ η OH . In the middle of the gap (k → 0), one retrieves the constant value pointed out by Bhowal and Vignale [16],

σ OH = e 2π m e v 2 F 3ℏ 2 g L 1 |∆ + λSOC | + 1 |∆ -λSOC | .
(13) This expression indicates that the orbital Hall conductivity decreases when increasing the gap, independently of the topological nature of the gap. To assess the interplay between the spin-orbit coupling strength λ SOC and the staggered potential ∆, Fig. 2(a) [Fig. 2(b)] shows the orbital Hall conductivity as a function of ∆ (λ SOC ) for different values of λ SOC (∆). We find that the orbital Hall conductivity reaches a maximum whenever ∆ ≈ λSOC , i.e., at the topological phase transition. We note that the orbital Hall conductivity tends to be larger in the trivial phase than in the topological phase, as confirmed by the phase diagram presented in Fig. 2(c). In this panel, brighter regions correspond to larger values of the orbital Hall conductivity, located in the topologically-trivial regions, |∆| > | λSOC |.

Finally, in order to establish a connection between the Kane-Mele model and the realistic two-dimensional materials discussed in the next section, we report the orbital Hall conductivity as a function of the gap size in Fig. 2(d). Our calculations confirm that the orbital Hall conductivity systematically decreases with the size for the gap, be it driven by ∆ or by λSOC consistently with Eq. ( 13). When the staggered potential is turned off (∆ = 0, red curve), at small λSOC , the gap is located close to the K and K' points, as shown in Fig. 3, and Eq. ( 13) applies. However, in the large spin-orbit coupling limit, i.e., λSOC ≈ t, the gap has moved to the M point (Fig. 3). In this case, increasing the spin-orbit coupling strength maintains a gapped spectrum while incrementing the Fermi velocity, which leads to an enhancement of the orbital Hall effect, as depicted in Fig. 2(d). The two situations reported in Fig. 3 are representative of the case of germanene (small λSOC , Dirac cones at K and K' points) and bismuthene (large λSOC , Dirac cone at Γ point) discussed below.

III. ORBITAL TRANSPORT IN REALISTIC 2D MATERIALS

We now turn to the simulations performed in realistic two-dimensional materials presenting different topological characters. For the trivial systems, we consider two cases: h-BN/graphene bilayer where the proximity effects lead to the breaking of inversion symmetry [34,35], as well as hydrogene-decorated graphene in which a colossal enhancement of the spin-orbit coupling has been predicted [36]. For the non-trivial systems, we have selected bismuthene, proven to be a topological insulator with a sizeable gap in its buckled hexagonal structure [37], and germanene, characterized by a buckled structure and a large enough spin-orbit coupling capable to open a topological gap [38,39]. For the DFT [40,41] simulations, we used the Perdew-Burke-Ernzerhof [42,43] exchange-correlation functional. We performed the relaxation with the plane-wave basis as implemented in the Vienna Ab-initio Simulation Package (VASP) [44,45], and employ a plane-wave expansion cutoff of 400 eV along with a force criterion of 0.2×10 -2 eV/ Å with a (15 × 15 × 1) k-points sampling of the Brillouin zone. The ionic potentials were described using the projector augmented-wave (PAW) method [46]. Finally, the Hamiltonian matrix was obtained through the Wannier90 package. The tight-binding representation was obtained by using a set of localized Wannier functions via the Wan-nier90/VASP interface [47]. To do so, we first provide a trial set of functions which represent the actual atomic orbitals in the system under study leading to a Hamiltonian written in terms of orthogonal Wannier functions. The evaluation of any physical observable is performed from the obtention of the eigenenergies and eigenvectors of this Hamiltonian [47].

h-BN/graphene and graphene+H

The recent proposal suggesting gapped graphene as an orbital Hall insulator has shed light on the nature of the valley Hall effect [16]. Motivated by this realization, we present the results on h-BN/graphene heterostructures. This system has been extensively studied in the last years in various contexts [34,48,[START_REF] Salvador-Sánchez | Generation and control of nonlocal chiral currents in graphene superlattices by orbital hall effect[END_REF]. While free-standing graphene is a topological semimetal that possesses a robust band structure protected by inversion symmetry, it loses this symmetry by proximity with a material like h-BN. The whole heterostructure resembles a Kane-Mele model in the trivial phase with a gap opening whose size is given by the interaction with h-BN. Theoretically it has been shown that spin manipulation would be possible in this scenario [35]. Most importantly for our purpose, graphene acquires a p x -p y orbital hybridization when interfaced with h-BN, which promotes the onset of intra-atomic orbital Hall effect. The Hall conductivities are shown in Fig. 4 where a small spin Hall effect is observed [Fig. 4(b)], whereas the orbital Hall response is one to two orders of magnitude larger. In particular, the intra-atomic orbital Hall effect [Fig. 4(c)] displays a moderate value within the energy window around the charge neutrality point while the total orbital Hall effect, which contains both intra-and inter-atomic contributions, [Fig. 4(d)] attains the largest value of the three Hall responses. Notice that the energy profile of the Hall responses are similar, as both spin and orbital Hall effects are driven by proximity with h-BN. We point out that d orbitals were previously introduced in order to increase the accuracy of the band structure compared with that obtained from GW+DFT simulations [START_REF] Boykin | [END_REF]. By enlarging the orbital basis to account for such d orbitals, we found that their contribution to the orbital Hall conductivity is close to ∼ 5%. This is rather negligible, especially in the present case where the largest contribution comes from the total orbital moment and is not related to the atomic orbital character of the Bloch state. We now consider graphene decorated with hydrogen. The inclusion of hydrogen is sufficient to enhance graphene's spin-orbit splitting up to 100 µeV locally [36,51]. For this system, we have considered a 5 × 5 supercell with a single hydrogen atom on top at the center of the graphene flake. The band structure is presented in Fig. 5(a) where spin-orbit coupling was also taken into account, showing a good agreement with previous reports [51]. The orbital and spin textures are shown in Fig. 5(b) and (c) respectively, for the most energetic valence bands closer to the hydrogen states, well localized in the spectrum. Our results suggest a large imprinted p x -p y hybridization which leads to the large value of the atomic orbital momentum L z having hot spots at inequivalent points in the hexagonal Brillouin zone. This is encouraging from the orbital transport perspective. Our calculations show that the intra-atomic orbital Hall conductivity (ACA) is around 0.6 (e/2π) whereas the total (intraand inter-atomic) orbital Hall conductivity is about 2.1 (e/2π), which is comparable to the h-BN/graphene case discussed above. In contrast, the spin Hall conductivity is about ≈ 0.15 (e/2π), still much smaller than the orbital Hall effect, but one order of magnitude larger than the spin Hall effect computed in h-BN/graphene, demonstrating the large spin-orbit coupling enhancement in this heterostructure. This result is remarkable especially considering that it solely arises from the interaction between graphene and hydrogen. On the other hand, the spin texture induces a local magnetic moment of ∼1 µ B , leading to a large spin splitting. The orbital Hall effect will also appear in other graphene-based heterostructures whose band structure is tuned by proximity effects [14,52,53]. The next material we consider is germanene which pos-sesses a narrow gap and has a buckled structure that favors a sp 3 hybridization inducing a p x -p y hybridization away from the neutrality point [Fig. 6(a)] which results in an orbital texture in momentum space [Fig. 6(b)]. Alike Kane-Mele model with small spin-orbit coupling, germanene possesses slightly gapped Dirac cones located at K and K' points [Fig. 6(c)]. The Hall conductivities are depicted in Fig. 6(d) where the spin Hall effect (blue) reaches a (narrow) quantized plateau at the Fermi level, associated with the non-trivial phase. Whereas the spin Hall conductivity is peaked close to the gap, where the spin Berry curvature is maximum, the non-vanishing orbital texture in germanene leads to a finite value of the orbital Hall conductivity (intra-atomic contribution in black, total contribution in green) on a much broader range of energy around the gap. Notice that the total orbital Hall effect remains smaller than the intra-atomic Hall effect, which implies that inter-atomic and intraatomic contributions partially cancel each other. A particular feature of germanene (and bismuthene, see below) is that inversion symmetry is preserved, and therefore the total orbital Hall response has its origin in the nonabelian nature of the Berry curvature as already shown [18]. The last system we consider is buckled bismuthene, which displays a much larger gap than germanene due to a much larger spin-orbit coupling [Fig. 7(a)]. In this crystalline phase, bismuthene's band character is inverted at Γ point due to spin-orbit coupling, following the same process as described in the Bernevig-Hughes-Zhang model [54,55], exemplified above by the Kane-Mele model with strong spin-orbit coupling (see Fig. 3). This change in the band structure allows for a strong s-character at this point in reciprocal space leading to a quenched orbital texture [56]. The spin, intra-atomic and total orbital Hall conductivities are displayed on Fig. 7(b). The absence of an orbital texture in terms of the p x and p y near the gap leads to a vanishingly smaller intra-atomic orbital Hall effect (black), which increases away from the gap due to enhanced p x -p y hybridization. In contrast, the total orbital Hall conductivity reaches a large value (green), even larger than that of the spin Hall conductivity (blue). These larger values can be understood by looking at the orbital moment distribution along the momentum path M -K -Γ -M shown in Fig. 7(c). The hot-spots located at the Γ-point lead to a larger value of the total orbital Hall conductivity compared to the spin one. Notice that the spin conductivity is quantized, whereas the total orbital conductivity is not. We mention in passing that it has been recently suggested that such non-quantized plateaus are related to highorder topological insulating behavior [START_REF] Costa | Connecting higher-order topology with the orbital hall effect in monolayers of transition metal dichalcogenides[END_REF]. Hence, the overall scenario in bismuthene contrasts markedly with that in germanene and follows the situation discussed in the previous section using the Kane-Mele model. Consequently, from a materials' perspective, we are able to draw differences on the orbital response based on details of their band structures and orbital character.

Germanene and bismuthene

IV. IMPACT OF DISORDER ON ORBITAL HALL TRANSPORT

An important question that remains unanswered at this point is the impact of disorder-induced scattering on the orbital conductivity. As a matter of fact, a simpleminded rationale suggests that the intra-atomic orbital Hall effect, which arises from the atomic orbital moment, would be less sensitive to momentum scattering than the inter-atomic orbital Hall effect, which arises from selfrotation of the electron wave packet in the unit cell. To investigate the impact of disorder, we consider three different systems: germanene, h-BN/graphene, the two narrow-gap semiconductors studied above, and MoS 2 , a large band gap semiconductor that has been predicted to be an orbital Hall insulator [6,[START_REF] Cysne | Rappoport[END_REF]. From the tightbinding basis obtained by ab initio simulations, we introduce disorder by the inclusion of an on-site Anderson disorder which can be expressed mathematically like

H = H 0 + i V i , (14)
where H 0 is the bare Hamiltonian corresponding to a 10 × 10 supercell and V i is an onsite potential acting on Hall conductivity over its maximum value σ 0 OH in the pristine case. We computed both the intra-atomic (ACA -red) and total contributions (green). The inset in (a) shows a sketch of a disorder realization in the lattice.

the i site with values [-1, 1]eV . We calculated the Hall conductivity for 40 random realizations of every fixed set of parameters.

The results obtained for germanene and h-BN/graphene are depicted in Fig. 8(a,b) and (c,d), respectively. To understand how disorder affects the orbital transport, we have considered two different transport regimes: (i) the single band case, where the carrier's energy to close to the gap, ε = 0.1 eV (solid lines), and (ii) the multiband case, where the carrier's energy is far from the gap, ε = -3.0 eV (dashed lines). In the former, the band dispersion is mostly linear and the longitudinal conductivity of germanene (a) and h-BN/graphene (c) slowly decays as a function of disorder due to enhanced scattering. In contrast, when the energy lies far from the gap, in the multiband case (dashed), the conductivity decay is more dramatic, as expected in conventional metals. We have also computed the intra-atomic and total orbital Hall conductivities for these different situations (b,d). To better visualize the effect of disorder, we report the ratio between the Hall conductivities with and without disorder, σ OH /σ 0 OH . In the single band transport regime, we find that the intra-atomic conductivity (red) is mostly flat, independent on the disorder. Nonetheless, the total orbital Hall effect (green), which contrains both intraand inter-atomic contributions, is as a whole much more sensitive to disorder and decreases continuously. In fact, the intra-atomic Hall effect is controlled by the orbital Berry curvature of the single band and is therefore expected to be rather robust against disorder whereas the inter-atomic Hall effect, which arises from self-rotation of the wave packet in the unit cell is much more sensitive to onsite energy fluctuations brought by Anderson disorder.

In the multiband transport regime (dashed lines), we find that both the intra-atomic and total orbital Hall effect decay at a similar rate. This distinct behavior suggests that the linear dispersion of the single band transport regime has a strong impact on the robustness of the intra-atomic orbital Hall effect. In contrast, the total orbital Hall effect, that contains the inter-atomic contribution, is much more sensitive to Anderson-type disorder. We must note that the intra-atomic contribution is larger than the total one in the case of germanene while the opposite is true for h-BN/graphene. This result indicates that the inter-atomic contribution is much more sensitive to disorder than the intra-atomic one, irrespective of the transport regime. We now turn our attention to MoS 2 , a large band gap orbital Hall insulator studied previously [6,19]. The band structure and orbital Hall conductivity are reported in insert for reference. Here, we fix the carrier's energy at ε=1.5 eV above the center of the gap in a region where the energy dispersion is quadratic. The value of the conductivity decreases rapidly with increasing the impurity concentration, as expected in a conventional, leading to a decay over about two orders of magnitude, quite different from linearly dispersing narrow-gap germanene and graphene/hBN. The behaviour of the (normalized) orbital Hall conductivities is depicted in Fig. 9(b) where dashed lines correspond to a transport energy of ε=1.5 eV and the solid lines correspond to a transport energy taken in the middle of the gap. The inset shows the absolute values for reference. The behavior we obtain is qualitatively similar to the one observed in the narrowgap semiconductors discussed previously. In the gap, the intra-atomic orbital Hall conductivity (solid red) is insensitive to disorder, as expected from a Berry-curvature driven effect, whereas the total orbital Hall conductivity decays. It has been argued recently that the in-gap intraatomic orbital Hall conductivity is associated with intraatomic orbital polarized edge states [START_REF] Cysne | Rappoport[END_REF]19] that remain insensitive to the disorder, although this picture might change when considering open boundary conditions like in nanoribbons for instance [START_REF] Pezo | [END_REF]. This contrasts to what we find for the total Hall response where a drop of nearly half of the initial value is observed. When the energy is set in the conduction band, one finds that both intraatomic and total orbital Hall conductivities decrease with a similar rate, as already observed in narrow-gap semiconductors.

To improve our study of the disorder we have considered 30 × 30 supercells of a two-dimensional model system representing d xy -d x 2 -y 2 and p z orbitals for a 3-band MoS 2 [59], capable to display both intra-atomic and inter-atomic orbital Hall effects (not shown). We were able to retrieve similar results calculated using the full Wannier Hamiltonian: in the gap, the total orbital Hall response decreases upon increasing disorder while the intra-atomic orbital Hall effect is preserved. These calculations lead to qualitatively the same conclusions as for the system with a larger number of orbitals. Our simulations therefore suggest that the total orbital Hall contribution is more sensitive to disorder and corroborates our previous comment regarding its relation with the periodic functions living in the bulk.

The present study, based on a real-space random potential distribution, in principle covers extrinsic scattering mechanisms such as side-jump and skew-scattering, which are known to be central to spin and anomalous Hall effects [60]. Indeed, the inclusion of disorder leads to a drop of the orbital Hall conductivity for states in the valence band, which seems to corroborate the conclusion drawn in [61]. In this work, it is shown that for a 2-band model system the side-jump contribution totally cancels the Berry-curvature contribution. Furthermore, our results exhibit a similar behavior for both the intra-atomic and the total orbital Hall contributions in the metallic regime.

V. CONCLUSION

In summary, we have explored the microscopic origin of the orbital Hall effect in model and realistic twodimensional Dirac materials. Since the orbital Hall effect is intimately connected with the Berry curvature of the material, we first investigated the inter-atomic orbital Hall contribution in the Kane-Mele model, that accommodates topological phase transition and in which the intra-atomic contribution is absent. We found that although the orbital moment itself behaves differently in the topologically trivial and non-trivial phases, the resulting orbital Hall conductivity is rather controlled by the size of the gap, irrespective of its topological nature.

We then studied the orbital Hall effect in selected twodimensional materials, starting with graphene. Whereas orbital Hall effect is absent of pristine graphene, it can be turned on by inducing a global or local gap, either interfacing graphene with h-BN or by using hydrogen adatoms, respectively. In these cases, the emergence of orbital and spin textures in reciprocal space stands out as a key ingredient for the generation of orbital and spin polarized Hall currents. These predictions are particularly intriguing given that these two systems are made out of light elements unable to portray a sizeable spinorbit coupling by their own. Although the experimental distinction between entangled spin and orbital signals remains a challenge even in materials with moderate spinorbit coupling, our results are encouraging as they show that light materials with negligible spin orbit coupling would display a large orbital Hall effect.

We then moved on to investigate the orbital Hall currents in two selected two-dimensional topological insulators, germanene and bismuthene, which represent two distinct realizations of the Kane-Mele model, with weak and strong spin-orbit coupling, respectively. In germanene, we have found that the intra-atomic orbital Hall contribution displays a larger value than of the total orbital Hall one, resembling the weakly spin-orbit coupled non-trivial phase of the Kane-Mele model (small spinorbit gaps at K and K' points). The existence of the quantum spin Hall effect in germanene is corroborated with a narrow plateau appearing for the spin Hall conductivity while the values for the orbital Hall effect re-mains larger in a broader energy window. In bismuthene, besides the spin Hall effect, we have found a large orbital Hall conductivity coming from the orbital moment carried out by the bands near the Fermi level. This situation resembles the strongly spin-orbit coupled non-trivial phase of the Kane-Mele model (large spin-orbit gap at M point).

Finally, we investigated the impact of disorder on the intra-atomic and inter-atomic contributions of the orbital Hall effect in two-dimensional systems featuring very different transport regimes (insulating, single-band and multiband metallic regimes). We find that the intraatomic orbital Hall effect tends to be less affected by disorder than the total orbital Hall effect, especially in the insulating and single-band regimes, i.e., in situations where the orbital Berry curvature is smooth and well-defined. In contrast, in the multiband transport regime, both intra-atomic contribution and total orbital conductivity. These results suggest that irrespective of the transport regime, the intra-atomic part of the orbital Hall effect is more robust than the inter-atomic part.

The present work sheds light on the mechanisms responsible for orbital Hall effect in two-dimensional materials, and in particular clarifies the role of the gap. The intimate connection between the orbital Hall transport and the Berry curvature of the band structure opens interesting perspectives for the external control of the orbital transport through interfacial engineering or strain, as demonstrated in Ref. [62]. From this standpoint Van der Waals heterostructures made of light materials, such as graphene and h-BN for instance, could be used for the realization of nonlocal orbital devices, akin to the all-electric valley or spin Hall transistor [63,64]. In this context, a more comprehensive understanding of the orbital relaxation induced by momentum scattering is necessary.

B

  FIG. 1. Kane-Mele model band structure for the trivial phase with the Berry curvature displayed in the color bar (a). The orbital moment plotted for trivial (∆ = 0.25t, λSOC = 0) (b) and topological (∆ = 0, λSOC = 0.25t) (c) phases for the two valence energy bands. The black (red) curve refers to the lowest (highest) valence band.

FIG. 2 .

 2 FIG. 2. (Color online) (a) Orbital Hall conductivity of the Kane-Mele model as function of ∆ for fixed values of the λSOC [0.25t (black), 1.0t (red) and 1.5t (green)]. (b) Orbital Hall conductivity as function of λSOC with fixed values of ∆ [0.25t (black), 1.0t (red) and 1.5t (green)]. (c) Phase diagram for the orbital Hall conductivity as function of both ∆ and λSOC. Whenever |3 √ 3λSOC| > |∆|, the orbital Hall effect is larger than for the opposite situation. (d) Orbital Hall conductivity in the center of the band gap as a function of ∆ (black) and λSOC (red).

FIG. 3 .

 3 FIG. 3. Band structure for the Kane-Mele model with different values of λSOC. In the non-trivial phase, the increasing of spin orbit coupling not necessarily leads to the increase of the gap, this is in connection to what is depicted in Fig. 2(d).

FIG. 4 .

 4 FIG. 4. (a) h-BN/graphene heterostructure electronic band structure showing the band gap opening at the Dirac point the inset, (b) spin Hall, (c) intra-atomic orbital Hall and (d) total orbital Hall conductivities. The grey shaded region in (a) corresponds to the energy window where the Hall conductivities were calculated.

FIG. 5 .

 5 FIG. 5. Graphene+H electronic band structure considering spin-orbit coupling (a). Orbital texture (b) and spin texture (b) for the most energetic valence band. In this case the isolated flat bands come from the Hydrogen atom.

FIG. 6 .

 6 FIG. 6. (a) Germanene projected density of states showing the pz (green), px , py (orange) and s (blue) states as a function of the Fermi level. (b) Corresponding orbital texture for the most energetic valence band. (c) Germanene electronic band structure, and (d) spin (blue), intra-atomic (black) and total orbital Hall conductivities (green).

FIG. 7 .

 7 FIG. 7. (a) Bismuthene electronic band structure. (b) Spin (blue), intra-atomic (black) and total orbital Hall conductivities (green). (c) Orbital moment calculated for the two most (red and black) and two less energetic (blue and green) bands along the M -K -Γ -M kpath.

FIG. 8 .

 8 FIG. 8. (a,c) Longitudinal conductivity and (b,d) orbital Hall conductivity near the Fermi level (ε = 0.1 eV) for germanene (a,b) and h-BN/graphene (c,d) as a function of the disorder concentration, the values plot correspond to the ratio orbital Hall conductivity over its maximum value σ 0OH in the pristine case. We computed both the intra-atomic (ACA -red) and total contributions (green). The inset in (a) shows a sketch of a disorder realization in the lattice.

FIG. 9 .

 9 FIG. 9. (a) Longitudinal conductivity for an energy of 1.5 eV with respect to the zero energy (middle of the gap) defined as the Fermi level. (b) Orbital Hall conductivities over the maximum value (σ 0 OH )in the pristine system, calculated by the Kubo formula in a 10×10 supercell containing Anderson type disorder, intra-atomic Hall conductivity (ACA) and total orbital Hall conductivity are depicted with dashed lines for an energy of 1.5 eV with respect to the zero energy defined as the Fermi level. The solid lines correspond to the Hall conductivites calculated the Fermi level. The inset shows the actual values of the orbital conductivity for both intra-atomic and total responses.

VI. ACKNOWLEDGMENTS

This work was supported by the ANR ORION project, grant ANR-20-CE30-0022-01 of the French Agence Nationale de la Recherche. D. G.O. and A. M. acknowledge support from the Excellence Initiative of Aix-Marseille Université -A*Midex, a French "Investissements d'Avenir" program.


Emilie Jué

C K Safeer

Marc Drouard

Alexandre Lopez

Paul Balint

Liliana Buda-Prejbeanu

Olivier Boulle

Stephane Auffret

Alain Schuhl

Aurelien Manchon

Ioan Mihai Miron
email: mihai.miron@cea.fr

Gilles Gaudin

Chiral damping of magnetic domain walls

Structural symmetry breaking in magnetic materials is responsible for a variety of outstanding physical phenomena. Examples range from the existence of multiferroics 1 , to current induced spin orbit torques (SOT) 2-7 and the formation of topological magnetic structures 8-11 . In this letter we bring into light a novel effect of the structural inversion asymmetry (SIA): a chiral damping mechanism. This phenomenon is evidenced by measuring the field driven domain wall (DW) motion in perpendicularly magnetized asymmetric Pt/Co/Pt trilayers. The difficulty in evidencing the chiral damping is that the ensuing DW dynamics exhibit identical spatial symmetry to those expected from the Dzyaloshinskii-Moriya interaction (DMI) 12-18 . Despite this fundamental resemblance, the two scenarios are differentiated by their time reversal properties: while DMI is a conservative effect that can be modeled by an effective field, the chiral damping is purely dissipative and has no influence on the equilibrium magnetic texture. When the DW motion is modulated by an in-plane magnetic field, it reveals the structure of the internal fields experienced by the DWs, allowing to distinguish the physical mechanism. The observation of the chiral damping, not only enriches the spectrum of physical phenomena engendered by the SIA, but since it can coexists with DMI it is essential for conceiving DW and skyrmion devices 19 .

The role of the SIA on the magnetic properties of the metallic ferromagnets has been the focus of extensive research. Whether the SIA occurs at an interface or in the bulk of a material, it may induce non-collinear exchange interactions 8 (DMI). This gives rise to a rich variety of magnetic phases, such as spin helices 9,10 and skyrmions 11 or to magnetic frustration 23 .

 this is an important idea, as it may solve a long-standing issue of DWs moving in the direction of the electric current and not into that of the electron flow, as one would expect

from spin transfer torque [START_REF] Thiaville | Dynamics of Dzyaloshinskii domain walls in ultrathin magnetic films[END_REF][START_REF] Ryu | Chiral spin torque at magnetic domain walls[END_REF][START_REF] Emori | Current-driven dynamics of chiral ferromagnetic domain walls[END_REF] . Moreover, the observation that the field-driven DW motion in Pt/Co based multilayers is made asymmetric by in-plane magnetic fields was also linked to the possible occurrence of DMI [START_REF] Je | Asymmetric Magnetic Domain-Wall Motion by the Dzyaloshinskii-Moriya Interaction[END_REF][START_REF] Hrabec | DMI meter: Measuring the Dzyaloshinskii-Moriya interaction inversion in Pt/Co/Ir/Pt multilayers[END_REF] .

Until now, DMI was the only well established interaction known to link the DW magnetic texture to the SIA of the layer structure. In this letter we propose a new mechanism, where the magnetic texture can sense the SIA as a damping whose intensity depends on the DW chirality.

DWs in materials with perpendicular magnetic anisotropy (PMA) are the perfect testground for such phenomena: they are chiral objects and their motion is relatively easy to detect.

From the perspective of the magnetization dynamics, all relevant physical phenomena translate into torques that control the magnetization according to the phenomenological Landau-Lifshitz-Gilbert (LLG) equation [START_REF] Gilbert | A phenomenological theory of damping in ferromagnetic materials[END_REF] . Their time reversal properties separate them in two categories. They can, either derive from a form of energy and be invariant to time reversal (anisotropy, exchange), or they may be dissipative and break time reversal symmetry (damping). These two scenarios are exclusive, in the sense that there is no other possibility.

Similarly, the physical phenomenon linking the DW motion asymmetry to the material's SIA must fall into one of these categories. It can only derive from chiral energy (DMI) or be a form of chiral dissipation.

The effect of chiral energy has already been evidenced experimentally [START_REF] Je | Asymmetric Magnetic Domain-Wall Motion by the Dzyaloshinskii-Moriya Interaction[END_REF] . The DMI implies the existence of an internal field HDMI setting the DW chirality (Figure 1b). As a consequence, the core magnetization of up/down (↑↓) and down/up (↓↑) DWs will point in opposite directions. The field driven DW velocity can sense this interaction if an external inplane field (Hip) is applied simultaneously with the perpendicular field (Hz). ↑↓ and ↓↑ DWs moving along Hip will experience different effective fields (Hip-HDMI and Hip+HDMI) and therefore will have different velocity. As long as the applied field does not alter the equilibrium DW profile (by significantly tilting the domains) the DW dynamics are correctly described by a collective coordinate model. Since in this case Hip and HDMI enter the DW equation of motion on equal footing [START_REF] Je | Asymmetric Magnetic Domain-Wall Motion by the Dzyaloshinskii-Moriya Interaction[END_REF] , DMI will bias the velocity dependence on Hip by +HDMI and -HDMI respectively (Figure 1b), making the motion asymmetric.

The chiral damping mechanism, although allowed by symmetry has never been observed experimentally. Its chief attribute is that, unlike DMI, it causes asymmetric displacements without biasing the magnetic configuration (Figure 1d).

We use wide field Kerr microscopy to probe the DW motion in asymmetric Pt30Å/Co6Å/PtxÅ trilayers subjected simultaneously to Hz and Hip. A pulse of Hz nucleates a reversed domain. A second pulse of Hz accompanied by Hip further propagates the circular DW. Figure 2a shows typical magnetic differential images corresponding to the DW displacement produced by Hz and Hip. The reversed magnetization produces bright (or dark) contrast. Without Hip the DWs move over equal distances in all directions. The direction and the magnitude of the asymmetry are controlled by Hip. This is qualitatively consistent with the existence of the HDMI: depending on the type of DW, (↑↓) or (↓↑), Hip either adds-up or subtracts to HDMI, leading to different effective field for the two DWs. To verify this possibility, we measured the DW velocity (vDW) vs. Hip for a fixed Hz. While Hip can be applied continuously, as it provokes no DW motion, Hz is applied in pulses of 400 ms. vDW is extracted by dividing the displacement observed on the Kerr image by the pulse length.

Figure 2 shows the effect of Hip on vDW. Besides the chiral effects emerging from SIA, Hip may also affect vDW through phenomena unrelated to SIA (such as the variations of the effective anisotropy). To separate the chiral effects we decompose vDW into a symmetric (S) and anti-symmetric (A) component with respect to Hip:

) 2 1 ( A S v + = ↑↓ and ) 2 1 ( A S v - = ↓↑ , where ( ) 2 ↓↑ ↑↓ - = v v S and ( ) ( ) ↓↑ ↑↓ ↓↑ ↑↓ + - ⋅ = v v v v A 2
. While S can be affected by non-SIA phenomena, A is a robust indicator of chiral effects, as by symmetry, only SIA related phenomena can create it. Since A dictates the ratio

A A v v - + = ↓↑ ↑↓ 2 2
, it is also an indicator of the asymmetry of DW motion the along the direction of Hip (Figure 2a). The variation of A vs. Hip has two distinctive features: it varies continuously up to a saturation field (Hsat=40mT), and then saturates (Figure 2a,2b), meaning that below Hsat the velocities of the two DWs vary differently, while above Hsat they vary identically to each other. This indicates that the parameter responsible for ↓↑ ↑↓ ≠ v v must also saturate at Hsat. In our samples, the only micromagnetic parameter that may vary monotonically up to 40 mT and then saturate, is the DW core magnetization. Due to the large uniaxial anisotropy (Hk=500 mT) 25 , the domain tilting affecting the DW internal structure and the effective DW width, becomes significant at much higher fields.

Hsat must then designate the field value for which both DWs become magnetized along Hip, implying that HDMI is smaller than Hsat. If DMI imposes homochiral Néel DWs, Hip acts only on the magnetization of one of the two DWs, as the other one is already aligned (Figure 1a). Therefore, vDW should display a feature at Hsat for only one of the DWs (Figure 1c). On the contrary, we observe (Figure 2c,S5) that vDW vs. Hip changes slope close to Hsat for both DWs.

The fact that Hsat is physically relevant for both DWs excludes the existence of homochiral Néel DWs in our samples. It also sets a ceiling value for HDMI within the precision of Hsat determination (±10mT).

As explained earlier, if A is produced by a small DMI (that we cannot rule out completely), the velocity curves for the two DWs should have the same shape and only differ through a horizontal offset given by the value of the HDMI bias field (Figure 1b). By shifting the two curves in opposite directions by HDMI, they should coincide. Clearly, the measured vDW do not follow this scenario: it is impossible to overlap the curves by any amount of lateral shifting (Figure 2c).

Based on the variations of vDW induced by Hip at constant Hz, we established that our samples support Bloch DWs and that the asymmetry mechanism evidenced in samples with Néel DWs [START_REF] Je | Asymmetric Magnetic Domain-Wall Motion by the Dzyaloshinskii-Moriya Interaction[END_REF][START_REF] Hrabec | DMI meter: Measuring the Dzyaloshinskii-Moriya interaction inversion in Pt/Co/Ir/Pt multilayers[END_REF] does not apply here. Now we switch perspectives and study the DW motion as a function of Hz at fixed Hip. To separate the effect of Hip on the chiral energy and chiral damping, we rely on the well-known creep scaling law, describing the DW velocity at very small driving fields:

        ⋅ - ⋅ = -4 / 1 4 / 1 0 exp z B p c creep H T k H U v v ( 1
)
c U -height of the pinning barrier; p H -pinning field;

Since the prefactor 31

      ⋅ ⋅ = T k U C f d v B c exp 0 0 0 ( 1 ≈ C
is an empirical constant) also depends exponentially on the pinning barrier, it is more convenient to write the creep scaling law [START_REF] Je | Asymmetric Magnetic Domain-Wall Motion by the Dzyaloshinskii-Moriya Interaction[END_REF][START_REF] Ferré | Universal magnetic domain wall dynamics in the presence of weak disorder C[END_REF] as

( )       - ⋅ - ⋅ = - - 4 / 1 4 / 1 0 0 exp c z B DW H H T k E f d v (2)
and group all the terms related to energy into a single exponential (

4 / 1 p c H U E = ). 4 C H H p c =
is the critical field that determines the limit of the creep regime 31 . Since the dissipation does not contribute to the exponent, the damping can only affect the prefactor 25-26 0 0 f d 31 (d0 -disorder correlation length; f0 -the attempt frequency).

As long as the creep scaling law applies, the plot of According to (2) they can be written ( )

        - - -         = ↓↑ ↑↓ - - ↓↑ ↑↓ T k E E H H f d f d A B mT

(

) ( ) ).

        - + - -         ⋅ = ↓↑ ↑↓ - - ↓↑ ↑↓ T k E E E H H f d f d f d S
Though it has been predicted 27 that friction affects the creep of an elastic membrane via the velocity prefactor, to our knowledge, in the case of DWs no microscopic mechanism has yet been proposed. For this particular case, the physical parameter most likely to depend on damping is the attempt frequency f0. We verify this possibility using micromagnetic simulations by calculating the DW depinning time from a single defect for different damping values (Figure 3c,d). Indeed, we observe a strong dependence of the release time on damping. Since the barrier height is kept constant in the simulations, it is only f0 that can be responsible for this variation.

We find that f0 is proportional to the inverse of the damping (DW mobility). This trend (

α 1 = DW v
) is further confirmed for a wider range of applied fields and damping values using a 1D model 25 .

The above mechanism can account for the chiral DW dynamics if the damping coefficient depends on the orientation of the DW's core with respect to the magnetization gradient:

( )

Z ip C m m ∇ ⋅ ∝ r α α
. However, this kind of damping cannot exist by itself in a real physical system since the dissipation may become negative if the chiral contribution exceeds the intrinsic damping. To prevent this unphysical behavior, the chiral damping must include at least a second component, which always offsets its value in the positive range. For example, its mathematical expression could be ( )

Z ip C m m ∇ ⋅ + ∝ r α α α 0 , where C α α > 0 .
Based on this expression, the amplitude of the chiral damping contribution can be extracted from the value of A. Above Hsat, where mip saturates,

c c v v α α α α - + = ↓↑ ↑↓ 0 0 and 1 1 2 / + - = ↓↑ ↑↓ ↓↑ ↑↓ v v v v A we obtain 0 2 / α α c A =
For the sample exhibiting the largest asymmetry (Figure S1), the amplitude of the chiral damping reaches approximately 50% of the constant part.

To illustrate the effect of the chiral damping on the DW dynamics we use a numerical collective coordinate model. Since such models do not include explicitly the magnetization gradients, we write the chiral damping using a simplified form:

φ α α α cos 0 C + =
. Here the value of the magnetization gradient is considered constant and is implicitly included in the values of α0 and αC. The orientation of the DW magnetization is described by the azimuthal angle φ. Figure 2e shows the computed vDW as a function of Hip at 0K and without pinning. The salient feature is that the velocities of both DWs saturate simultaneously, at the field value where their core magnetization saturates. In order to compare the experimental results with the simulations, the symmetric component S of the measured velocity must be removed, as it is not produced by the chiral damping. For this purpose we normalize v↓↑ and v↑↓ to the average velocity along x and y, the in-plane directions. The normalized DW velocity (Figure 2d, Figure S5) saturates at Hsat, in perfect agreement with the simulations.

At this stage, we do not understand why, contrary to previous studies, the effect of the chiral damping in our samples is larger than the effect of DMI. We infer that these two phenomena should co-exist, but their relative strength may be strongly dependent on the precise crystallographic structure. An independent hint of the co-existence of chiral damping and DMI was provided by the magnon dynamics in Fe/W bi-layers, where it was observed that the not only the dispersion relation was shifted due to DMI, but also that magnon lifetime depends on the chirality [START_REF] Zakeri | Magnon lifetimes on the Fe (110) surface: The role of spin-orbit coupling[END_REF] .

While theoretical work is definitely required to understand the origin of the chiral damping, we would like to speculate and point towards a plausible explanation: both theoretically and experimentally it was shown that the spin orbit torques have two components: a field like 2,3 (conservative) part and a damping like (dissipative) part [START_REF] Miron | Perpendicular switching of a single ferromagnetic layer induced by inplane current injection[END_REF][START_REF] Liu | Spin-torque switching with the giant spin Hall effect of tantalum[END_REF][START_REF] Fang | Spin-orbit-driven ferromagnetic resonance[END_REF] . Theoretically the DMI was linked to the conservative (field like) part [START_REF] Freimuth | Berry phase theory of Dzyaloshinskii-Moriya interaction and spin-orbit torques[END_REF][START_REF] Kim | Chirality from interfacial spin-orbit coupling effects in magnetic bilayers[END_REF] , as they share a common microscopic origin. In a similar way, the chiral damping reported here, might be reminiscent of its dissipative component.

From a practical perspective, the DW motion asymmetry can be useful for memory [START_REF] Parkin | Magnetic domain-wall racetrack memory[END_REF] or logic devices 32 requiring unidirectional DW movement. The application of an alternating Hz causes the domain to enlarge and then shrink to its initial position. When Hip is applied along with Hz, the DW motion will become asymmetric. If Hip changes sign together with Hz, the asymmetry of the domain growth will be opposite to that of its contraction. Consequently, the domain will not come back to its initial position but will be shifted along the axis of the inplane field [START_REF] Je | Asymmetric Magnetic Domain-Wall Motion by the Dzyaloshinskii-Moriya Interaction[END_REF] . The relative sign of the two field components controls the direction. Compared to current induced DW motion where trains of DWs are pushed along 1D wires, the phenomenon that we describe affords a distinctive ability of shifting 2D magnetic domains in any directions of a full sheet magnetic layer (Figure 4c; supplementary movie).

In conclusion, the quantitative analysis of field induced DW motion has revealed the existence of a chiral damping term whose effects on the DW dynamics have the same spatial symmetry as those expected from the DMI. Understanding and controlling the chiral damping is essential for the designing spintronics devices that rely on chiral magnetic textures, such as DW or skyrmion racetracks. b Without SIA and hence with no DMI to differentiate the two DWs, Hip will have the same effect on the field driven velocity of the two DWs that move along x (black dashed lines). If present, the HDMI bias field will shift laterally the velocity curves for the two DWs (red and blue lines). An important consequence is that any feature of the initial velocity curve allowing to identify the DW transformation, such as a change of slope when its magnetization saturates, will translate laterally, and always manifest at the same velocity (vsat). c. If no DMI is present, the DWs have Bloch structure; the structural changes produced by Hip will be identical for the two DWs. d. The chiral damping does not shift laterally the initial "non-SIA" velocity curve (black dotted line), but creates an asymmetric vertical correction. In this case, the DW transformation from Bloch to Néel, will not occur at the same velocity, but at the same field value (Hsat). nm wide nanowire pinned to a localized defect; the dark rectangle designates an area with 50% smaller anisotropy that attracts the DW d. The simulated depinning probability (60 independent events) for two damping values (0.5 and 0.25) at Hz=18mT. The de-pinning field at 0K is 30mT.

The inset shows that the curves overlap when the timescale is normalized by the damping constant. This indicates that f0 is proportional to the inverse of the damping. The variation of the anisotropy field reflects a change of the interfacial anisotropy. While Hk varies by 30% (0.55 T to 0.38 T) after correcting for the layer's demagnetizing field (µ0Ms=1.3 T), we calculate a maximum variation of the anisotropy constant of 10%.

Though the relative change of the interfacial anisotropy is small, it is nevertheless surprising to observe such a variation for Pt capping layers exceeding several nm. A plausible hypothesis is that the heating of the sample during deposition may influence its crystallographic structure in the vicinity to the interface. It is well known that high temperature annealing enhances the Pt and Co intermixing. In order to corroborate this hypothesis we have grown several series of the following samples: Pt30Co6Pt20, Pt30Co6Pt60 and Pt30Co6Pt3520. While the Pt60 layer of the second was deposited continuously, the Pt capping of the third was deposited in three steps consisting of Pt20 separated by 30 seconds pauses, to allow cooling. N.B. The changes of Hk correspond to small variations of the interfacial anisotropy constant, and therefore are sensitive to small variations of the growth conditions.

Despite the variations from series to series, we consistently find that Hk(Pt30Co6Pt20) > Hk(Pt30Co6Pt3520)+50mT > Hk(Pt30Co6Pt60)+100mT, in agreement with the different heating that the samples might have suffered.

To verify experimentally that the heat-induced changes in the anisotropy do not affect the asymmetry we annealed the tPt=18.7Å sample at 200°C in vacuum for 30 minutes, and measured again its Hk as well as the DW motion asymmetry. Indeed, we observe that the annealed sample has a smaller anisotropy field (the same as the sample with tPt=30Å), consistent with the Pt-Co inter-diffusing, while the asymmetry remains almost unchanged (Figure S1 d).

The only consequence of the anisotropy variation is that the DW width might vary from sample to sample. However, since it depends on the square root of the anisotropy field, its variations will not exceed 10%. Consequently, the DW demagnetizing field stabilizing the Bloch DW structure will vary by the same ratio. In the absence of DMI, Hsat is given by the DW demagnetizing field. The fact that we do not detect any significant sample-to-sample variation of Hsat is consistent with this scenario.

S2. The prefactor of the DW creep scaling law

Recently, Gorchon et al. [START_REF] Manchon | Theory of nonequilibrium intrinsic spin torque in a single nanomagnet[END_REF] evidenced experimentally an exponential dependence of the creep velocity prefactor on temperature. They found that:

      ⋅ ⋅ = T k U C f d v B c exp 0 0 0 0 0 f d
, is the product of the pinning potential correlation length (d0) and the depinning attempt frequency (f0). ( )

        - ⋅ - ⋅ = - - 4 / 1 4 / 1 4 / 1 0 0 exp c z B p c DW H H T k H U f d v
, where have two effects. On one hand, the slope is modified, but the curve is also shifted laterally, due to the variation of the pinning field (related to the energy barrier). On the contrary, the attempt frequency does not produce any variation of the slope; it only shifts the curve vertically. Our experimental data (Figure 3 in the main text) does not exhibit any significant change of slope when Hip is reversed (from -80mT to 80mT); we only observe a shift. We deduce that this must be a vertical shift (associated to f0), since a lateral shift (associated to the pinning field) would be accompanied by a change of slope. The counter example can be seen when comparing the curve at Hip=0mT to the curves obtained at 80mT. Their difference in slope is associated to a change of the pinning field (Figure 3 in the main text).

Order of magnitude of the attempt frequency.

Besides the measurements of Gorchon et al. [START_REF] Manchon | Theory of nonequilibrium intrinsic spin torque in a single nanomagnet[END_REF] the order of magnitude of d0 f0 can be extracted also from previous measurements of DW creep in Pt/Co/Pt samples with varying Co thickness [START_REF] Fang | Spin-orbit-driven ferromagnetic resonance[END_REF] . In this case it is the change of thickness that leads to changes of the effective parameters (exchange, anisotropy, pinning…) affecting the energy landscape. Figure 2d By considering that the disorder correlation length should be larger than the typical grain size (20 nm), the order of magnitude of the attempt frequency could be anywhere between 10MHz and 1GHz.

S3. Influence of damping on the DW velocity in the thermally activated regime

Thermally activated DW motion consists of relatively fast DW displacements interrupted by pinning and de-pinning events. If the driving field is sufficiently small compared to the pinning field, the DW spends most of the time trapped at pinning sites, trying to de-pin.

When the temperature fluctuations finally free it, the DW moves to the next pinning center. In this motion regime the experimentally determined DW velocity is much smaller than the steady DW velocity (which is effective in-between pinning events); it is merely a measure of the average de-pinning time. For this reason, instead of trying to assess the complete influence of damping on the thermally activated DW velocity we will first focus on the DW depinning.

The thermally activated de-pinning of a DW depends on two ingredients: the usually dominant one is the exponential dependence on the height of the energy barrier in units of thermal energy. The other, often neglected, is the pre-factor of this exponential, given by the attempt frequency ( = • exp ∆ / ). While the energy barrier does not depend on the damping value, the attempt frequency may be altered by variations of the damping. We analyze this possibility numerically.

a) Chiral damping and thermal fluctuations

The effect of temperature along with the chiral damping proves to be impossible to include in the numerical modeling. We find that the temperature alone induces DW displacements. This is a direct consequence of a combination of two factors:

1) the random thermal field is created using a pseudo random number generator;

2) the chiral damping makes the DW behave as a ratchet.

A ratchet is an object having an asymmetric response to excitations [START_REF] Miron | Perpendicular switching of a single ferromagnetic layer induced by inplane current injection[END_REF][START_REF] Liu | Spin-torque switching with the giant spin Hall effect of tantalum[END_REF] . When subjected to periodic excitations a ratchet produces a net directional motion (rotation, displacement, etc).

In the case of DW motion, the asymmetric response is given by the chiral damping. And since the temperature fluctuations are generated using a pseudo random number generator whose autocorrelation is not null, it includes periodic components that inevitably create an artificial displacement.

Despite this difficulty, it is still possible to test the effect of damping in the thermally activated DW motion, using an alternative approach. As long as mip is established by the competition between Hip and the DW's demagnetizing field, the only variable determining the damping value is Hip. Therefore, instead of computing the velocity dependence on Hip we can use the damping value as the free parameter. We use a constant damping (independent of mip)

and vary its value. These two methods are equivalent as long as Hz is sufficiently small, not to induce any supplementary DW deformation and mip is determined only by the balance between Hip and the DW's demagnetizing field. This approximation is valid in the case of very slow thermally activated DW motion (corresponding to our experimental conditions).

b) Depinning time vs. Gilbert damping

In the micromagnetic simulations (Figure 3d of the main text) we observed a net dependence of the depinning time on the damping coefficient. Since the energy barrier for the two simulations is the same, this difference may only originate in the DW's attempt rate.

Furthermore, the fact that the depinning time is proportional to the damping value shows that the DW mobility plays a significant role in the motion of the DW inside the pinning center. The DW depinning is a probabilistic process whose full characterization requires averaging over a statistically significant amount of events. Consequently, micromagnetic modeling, being relatively time consuming, is not the best-suited tool for the extensive study of the thermally activated DW motion. In order to further analyze the influence of the damping coefficient on the DW depinning we use a numerical collective coordinate model.

c) Velocity dependence on damping in the thermally activated regime

While 2D micromagnetics allow including realistically the partial pinning of the DW, by definition in the 1D models the entire DW is enclosed in the pinning potential. As a consequence, when projecting the 2D DW dynamics onto a single dimension, the effective pinning potential needs to be adapted accordingly. The simplest way is to use effective pinning centers that are much wider than the DW, such that the DW will have to shift far from the energy minimum before actually applying pressure on the pinning potential. The velocity calculated for a potential half the width of the DW (5nm). In this case, the DW velocity is independent of the damping value.

In order to probe the influence of the damping variation on the thermally activated DW motion, we calculate the velocity (Figure S3) as a function of the perpendicular magnetic field,

for several damping values. The DW parameters used in the simulations are ∆DW=5nm,

Hdip=40mT. The effective pinning field is modeled by a sinusoidal function with a period of 50nm, (10 times larger than the DW width) with amplitude of 3 mT. This form of pinning is chosen only for ease of implementing in the numerical model. Though it is not realistic for the quantitative calculation of the DW velocity, since it allows modeling the repeated pinning and de-pinning of the DW, it is sufficient to evidence the impact of damping on the DW velocity.

We observe that the DW velocity depends strongly on the damping value (0.2, 0.4, 0.6).

Moreover, by normalizing all the velocity curves by the inverse of the damping (proportional to the DW mobility) they overlap (figure S2a). This indicates that the DW velocity depends on the intrinsic DW mobility. On the contrary, when trapping the DW in a narrow potential well with the same barrier height (period of 2.5 nm and amplitude of pinning field of 60 mT) the velocity curves completely loose the damping dependence (figure S2b).

d) Attempt frequency in the 1D model

To identify the physical parameter responsible for the damping dependence of the DW velocity, we analyze the motion of the DW driven only by temperature (without Hz) inside the To verify this prediction, we repeated the measurements for different values of Hz.

Indeed we find that, though the velocity changes by three orders of magnitude, both Hsat and A remain the same. On the contrary, the symmetric component varies significantly, changing the overall shape of the curves. The uncorrelated behavior of the symmetric and asymmetric components of the velocity points to their different origin: while the asymmetry is given by the chiral damping through the attempt frequency, the symmetric part is linked to changes in the pinning potential induced by Hip.

S5. Link between the asymmetry and mip variation

The DW velocity measurements correspond to the creep regime, where the DW velocity depends linearly on the attempt frequency and exponentially on the DW energy and pinning strength. For simplicity, we merge all energy contributions under a same parameter (E)

( )

E f v exp 0 ⋅ ≈
We can write explicitly for the two DWs:

∂ ∂θ

The three parameters that define this energy and that can be in principle affected by Hip are: the magnetization gradient, the DW width and the mip orientation. Since the variation and the saturation of the asymmetry occur at low fields, while the magnetization gradient varies significantly at much larger fields, comparable with the uniaxial anisotropy, we deduce that the chiral energy variation is essentially due to a variation of mip. Therefore, independently of the mechanism that creates the asymmetry, its saturation at low fields, is a robust indicator of the mip saturation.

  of Hip on the prefactor and the exponent, as they appear as the intercept and the slope of the linear dependence (FigureS2).

  Figure 3a,b shows typical results for two samples with different capping thickness (Pt15.6Å and Pt30Å). The linear dependence specific for the creep regime is observed for the two samples and both DW polarities. To evidence the symmetric and anti-symmetric effects of Hip we consider

  Figure3a,b), implying that the bubble asymmetry is not due to chiral energy ( of non-chiral energy contributions related to Hip (

Figure 1 .

 1 Figure 1. Graphic illustration of the asymmetric DW dynamics in PMA materials. a. A DW

Figure 2 ..

 2 Figure 2. DW velocity modulation by Hip in Pt/Co/Pt layers. a. Kerr microscopy images of

Figure 3 .

 3 Figure 3. Effect of chiral damping on thermally activated DW motion a,b. Creep scaling law

Figure 4 .

 4 Figure 4. Bubble shifting by magnetic fields. Sequence of direct Kerr images showing the

Figure

  Figure S1. a. Asymmetry curves extracted from Kerr microscopy images for Pt thickness ranging from 15 Å to 70 Å. b. When normalized to 1 all the curves overlap, indicating that Hsat is the same for all the samples. The inset shows the dependence of the normalizing factor (saturation value of the asymmetry) with the thickness of the Pt capping. c. Measurement of the anisotropy fields of all the samples used in the experiment. d. Asymmetry of the annealed sample compared to the as deposited ones. If the asymmetry value would be correlated to the anisotropy field the asymmetry of the annealed sample (red) would be comprised in-between the asymmetries measured for the samples with 30Å and 56Å of Pt capping.

  constant. Inserting this formula into the creep scaling law leads to:

  By rewriting the creep law under this form all the energy contributions are removed from the velocity prefactor; they are included in Hc -1/4 , which is an offset of the abscissa.

Figure

  Figure S2 represents schematically the changes of the creep scaling law produced by variations of the energy barrier and of the attempt frequency. Changes of the pinning barrier

Figure S2 .

 S2 Figure S2. Influence of energy barrier and attempt frequency on the DW creep scaling law. a. The variations of the pinning barrier lead to a change of slope (red dotted line), but also to a lateral shift due to a change of the pinning field. The horizontal gray line represents the approximate transition from creep to flow. The black square is the approximate Hc -1/4 b. The variation of the attempt frequency leads only to a vertical shift of the curve. c. Hc -1/4 and d0 f0 can be roughly estimated from the intercept of the creep velocity curves corresponding to different pinning barriers.

  of supplementary reference 6, exhibiting the DW velocity in the 4 that all the curves corresponding to the different samples converge close to a small area. This area allows estimating the range of Hc -1/4 and d0 f0 values for the series of samples. Since our samples are nominally identical, and are made under the same conditions in the same sputtering chamber as these previous studies[START_REF] Fang | Spin-orbit-driven ferromagnetic resonance[END_REF][START_REF] Manchon | Theory of nonequilibrium intrinsic spin torque in a single nanomagnet[END_REF] , we use their estimate for the order of magnitude of estimate that the deviations from the creep regime should occur when 5

Figure S3 .

 S3 Figure S3. DW velocity in the 1D model for 3 different damping values. Each velocity point is obtained after averaging 30 displacements produced by 500ns long pulses a. DW velocity when the width of the pinning potential is 10 times larger than the size of the DW (5nm). The inset shows the velocity curves normalized by the inverse of the damping value. b. The velocity calculated for a potential half the width of the DW (5nm). In this

Figure S5 .

 S5 Figure S5. Asymmetry dependence on Hz. a. Velocity curves for different values of Hz. b. The corresponding asymmetries Araw = 2(v↓↑ -v↑↓)/(v↓↑ + v↑↓). c. The normalized velocities.

  below Hsat the two DW velocities vary differently, while above, their ratio the attempt frequency and energy variations coincidently compensate each other. This scenario can be excluded based on the robustness of the observation: the asymmetry saturation is observed for a wide range of Hip for all our samples. There are two possibilities left: Case 1) Damping induced asymmetry This mechanism links directly the two DW velocities to mip. Therefore, as indicated in the manuscript the asymmetry saturation coincides with the saturation of the DWs magnetization.Case 2) energy induced asymmetryIn thin films with vertical SIA and isotropic in the plane, the chiral energy for the two DWs is:

  for up/down and down/up DWs, and cos(ϕ) variations induced by Hip have the same sign, any variation of mip for any one of the DWs will amplify if mip is constant for both DWs.

† Currently at NIST Boulder, Colorado, US

Acknowledgements

We thank H. W. Lee, K-J Lee and T. Giamarchi for helpful discussions as well as A Thiaville, S Pizzini, and J Vogel for critically reading the manuscript and discussing the results. This work was partially supported by the ANR (11 BS10 008, ESPERADO) project and European Commission under the Seventh Framework Programme (GA 318144, SPOT) and (GA-2012-322369, Sport for Memory). AM has been supported by King Abdullah University of Science and Technology.

Competing financial interests

The authors declare no competing financial interests To probe the dependence of this phenomenon on the material properties, we have compared the asymmetries measured on several samples having different thickness of the Pt upper layer (15Å-70Å). We chose this approach, to ensure that neither the magnetic properties of the Co layer nor the interface asymmetry 1,2,3 , known to be essential for the DMI, are affected.

All the samples exhibit the same qualitative behavior: A varies up to Hsat, and fully saturates above. Surprisingly, A decreases as the top Pt is made thicker. The possibility to experimentally turn off A affords another independent way of probing its origin. The DMI creates an asymmetry by inducing a bias field that shifts the velocity curves towards positive Hip for one DW and towards negative Hip for the other (Figure 1c). Consequently, the saturation value of A should be proportional to the bias field (HDMI). And since Hsat is the experimental indicator of HDMI, the amplitude of A should be related to Hsat. Experimentally we do not observe such a connection: the variations of A are not associated to changes of Hsat (Figure S1), thus showing that A does not originate from HDMI.

The uniaxial anisotropy fields (Hk) were determined from measurements of anomalous Hall effect for all the samples. The results plotted in figure S1 shows a reduction of the anisotropy fields from the sample with the thinnest (15.6Å) Pt capping to the thickest (70Å). pinning potential. To estimate the DW attempt frequency we compute the frequency spectrum (by FFT analysis) of the DW position. We observe that the attempt frequency is not described by a well-defined peak, but rather by a cut-off value. When normalizing the different FFTs to the inverse of the damping (DW mobility) the spectra overlap. This proves that the characteristic attempt frequency of the DW (cut-off frequency) changes proportionally to its mobility. In the case where the pinning potential is not sufficiently wide (of the order of the DW width or smaller) the FFT shows a relatively welldefined resonance peak (corresponding to the DW's ferromagnetic resonance) whose position does not depend on damping. The fact that the attempt frequency (in this case given by the position of the peak) loses its damping dependence, explains why the DW velocity in the case of narrow pinning potential does not depend on damping.

S4. Asymmetry dependence on DW velocity.

The numerical analysis of the DW motion reveals an important feature: since the velocity variation induced by the damping does not depend of the value of Hz (Figure S2a), the measured asymmetry should also be independent of Hz.