Eric Cancès
email: eric.cances@enpc.fr

Paul Cazeaux
email: pcazeaux@umn.edu

Mitchell Luskin
email: luskin@umn.edu

Generalized Kubo Formulas for the Transport Properties of Incommensurate 2D Atomic Heterostructures

 to study aperiodic solids (disordered crystals, quasicrystals, and amorphous materials), notably in the presence of magnetic fields (quantum Hall e ect). We also present numerical approximations and test our methods on a one-dimensional incommensurate bilayer system.

Introduction

The synthesis and modeling of layered two-dimensional atomic heterostructures is currently being intensely investigated with the goal of designing materials with desired electronic and optical properties [START_REF] Geim | Van der Waals heterostructures[END_REF]. These multilayer two-dimensional materials are generally incommensurate, that is, the multilayer system does not have a periodic structure although each individual layer does have a periodic structure.

Despite the scientific and technological importance of incommensurate materials, an exact formulation has not yet been given for important properties such as the electronic density of states or the electrical conductivity. The typical approach to modeling and computing the properties of incommensurate structures is to approximate them by commensurate structures or supercells [START_REF] Terrones | Bilayers of transition metal dichalcogenides: Di erent stackings and heterostructures[END_REF]. Although this approach might provide a good approximation in many cases, the error is generally uncontrolled (see [START_REF] Cazeaux | Analysis of rippling in incommensurate one-dimensional coupled chains[END_REF] for an analysis of the supercell approximation of the mechanical relaxation of coupled incommensurate chains). Further, the approximation of the important case of small angle rotated bilayer structures requires supercell sizes that are too large for numerical solution [START_REF] Tritsaris | Perturbation theory for weakly coupled two-dimensional layers[END_REF]. In this paper, we give such an exact formulation and develop numerical approximations.

Within independent electron or mean field models such as Hartree-Fock or Kohn-Sham, the electronic density of states and transport properties for periodic structures can be rigorously formulated by the use of the Bloch transform to obtain generalized eigenstates of the Hamiltonian as a plane wave with wave vector in the Brillouin zone (k-space) multiplied by a periodic function. This approach leads to the classical Kubo formulae for the transport properties of periodic solids which can be formulated as the trace of corresponding operators [START_REF] Kaxiras | Atomic and Electronic Structure of Solids[END_REF].

Although incommensurate multilayer two-dimensional materials no longer have a periodic structure and the local environment of each atom is unique, this local environment can be simply characterized by shifts of each layer. Further, the shifts corresponding to the local environments are uniformly distributed over the periodic unit cell of each layer. One can then develop generalized Kubo formulae for incommensurate heterostructures by considering the integral over the uniformly distributed shifts of the trace of operators that depend of the corresponding local environment.

To give a precise formulation and enable mathematical and numerical analysis, we apply the C ˚-algebra approach for aperiodic solids introduced by Bellissard and collaborators [START_REF] Bellissard | Coherent and dissipative transport in aperiodic solids[END_REF][START_REF] Bellissard | The noncommutative geometry of the quantum Hall e ect[END_REF] to incommensurate heterostructures. Following these lines, we can characterize the incommensurate structure by its hull, which is a compact description of the local environments. We first present the construction of the hull for perfect multilayers in Sections 2.1 and 2.2. The case of disordered heterostructures is dealt with in Section 2. [START_REF] Bellissard | The noncommutative geometry of the quantum Hall e ect[END_REF].

In Section 3, we present the C ˚-algebra formalism for tight-binding models. Rather than working with operators on the infinite dimensional tight-binding state space, the C ˚-algebra approach allows us to exploit the simplification of working directly with a C ˚-algebra of functions which represent the operators of interest. Within the tightbinding model that we consider, we obtain a compact parametrization of the tight-binding Hamiltonian by using environment-dependent site and hopping functions [START_REF] Fang | Ab initio tight-binding hamiltonian for transition metal dichalcogenides[END_REF]. The C ˚algebra approach then allows us to concisely construct transport operators and their trace by utilizing the algebra structure. After introducing the abstract setting in Section 3.1, we describe the special case of perfect incommensurate bilayers in full detail in Section 3.2.

In Section 4, we present a new, minimalistic one-dimensional toy model to showcase the expected e ects of incommensurability in coupled multilayered systems. We chose to introduce here a simple discretization based on periodic supercells and on the Kernel Polynomial Method [START_REF] Weiße | The kernel polynomial method[END_REF]. In future works, we will develop more sophisticated approaches based on the C ˚-algebras introduced in this paper. We will present their numerical analysis, using the C ˚-algebra formalism, and use the minimalistic model introduced in this paper as a benchmark to analyze and test di erent implementation strategies, possibly targeting directly incommensurate cases, for computing the density of states, the conductivity or other observables. One such possible strategy is to exploit locality [START_REF] Carr | Twistronics: Manipulating the electronic properties of two-dimensional layered structures through their twist angle[END_REF][START_REF] Massatt | Electronic density of states for incommensurate layers[END_REF].

Previous research on the development of numerical methods to approximate transport properties within the C ˚-algebra formulation has been done by Prodan [START_REF] Prodan | Quantum transport in disordered systems under magnetic fields: A study based on operator algebras[END_REF] for the e ects of disorder and magnetic fields. Recent work on the analysis and computation of the density of states for incommensurate layers from the operator point of view is given in [START_REF] Massatt | Electronic density of states for incommensurate layers[END_REF]. We note that the density of states is defined in [START_REF] Massatt | Electronic density of states for incommensurate layers[END_REF] as a thermodynamic limit, while the density of states and transport properties are given explicit expressions in the C ˚-algebra approach, which directly provide the values of the quantities of interest in the thermodynamic limit.

Geometries of multi-layered systems

Perfect multilayer structure

Heterostructures are vertical stacks of a few two-dimensional crystalline monolayers. An example is depicted in Figure 1. Five distinct layers with di erent atomic components and structures are positioned as a vertical stack. Due to the weak van der Waals nature of the interactions between these layers, they do not relax into a common periodic structure, but rather each layer essentially keeps the structure it possesses as an isolated monolayer. The resulting assembly is thus in general not periodic: it is an aperiodic structure with a long-range order. A systematic model of heterostructures starts naturally by a rigorous depiction of this particular geometry, which can be idealized in the following way.

For the sake of generality, we consider d-dimensional systems embedded in pd `1qdimensional space. Note that we choose d as the natural dimension of the structure, since the pd `1q st dimension plays a very particular and limited role, especially in tightbinding models. We are in particular interested in d " 2 for the layered heterostructures which motivate our study, while the choice d " 1 enables us to present simple numerical examples in Section 4.

Let us consider such a pd `1q-dimensional system of p parallel d-dimensional periodic atomic layers denoted L j Ä R d`1 , j " 1 . . . p. We denote by pe 1 , ¨¨¨, e d`1 q an orthonormal basis of the physical space such that each layer is perpendicular to e d`1 ; from now on we identify the physical space with R d`1 using the cartesian coordinates x " px 1 , ¨¨¨, x d`1 q T associated with this basis; -h j the pd `1q st coordinate of the center of layer j. Without loss of generality, we can assume that 0 "

h 1 † h 2 † ¨¨¨ † h p ; -R j the d-dimensional periodic lattice of layer j;
-E j the matrix in R dˆd whose columns form a basis generating the lattice R j :

R j :" E j Z d Ä R d ; ( 1 )
j :" R d {R j the quotient of R d by the discrete lattice R j , which has the topology of a d-dimensional torus and can be canonically identified with the periodic unit cell p j :" E j r´1{2, 1{2q d of layer j; -m j the motif of layer j, i.e., the measure on R d`1 supported in p j ˆR and representing the nuclear distribution fl nuc j of layer j. More precisely, m j is a finite sum of positively weighted Dirac measures of the form

m j " M j ÿ k"1 z pjq k " x pjq k ,
where M j is the number of nuclei per unit cell in layer j, z pjq 1 , ¨¨¨, z pjq M j the atomic charges of these nuclei, and x pjq k P p j ˆR their positions in a reference configuration in which the center of layer j belongs to the plane x d`1 " 0. The nuclear distribution fl nuc j of layer j is then given by

fl nuc j " ÿ n j PR j m j `¨´p" j `nj `hj e d`1 q ˘" ÿ n j PR j M j ÿ k"1 z pjq k " " j `nj `xpjq k `hj e d`1 ,
where " j P p j depends on the horizontal position of the lattice sites of layer j relatively to the origin of the coordinates. Here and in the sequel, we use the same notation to denote a vector y of R d and its canonical embedding py T , 0q T in R d`1 . Remark 2.1. Note that, by virtue of the identification j " p j , an equivalence class " j P j can be seen either as a discrete set of points in R d , or as one point of the periodic unit cell p j . Since both viewpoints are employed here, we denote these two usages di erently to avoid confusion:

-" j will be used to denote a single point of the periodic unit cell p j , -" j `Rj will be used to denote the corresponding set of points in R d . In particular, the sum of the values of a function f over the lattice sites of layer j will be denoted as ∞

p j P" j `Rj f pp j q.

Translation group and the hull

To understand the geometry of our heterostructures, it is a useful exercise to picture the origin of coordinates as our viewpoint (in the sense of "position of observation").

Neighboring atoms then constitute a local environment, e.g., A-A stacking (aligned layers) vs. A-B stacking (staggered layers) in a graphene bilayer [START_REF] Fang | Electronic structure theory of weakly interacting bilayers[END_REF]. This environment, i.e., the positions of all atoms relative to the origin, constitutes a choice of configuration for the structure.

In the case of a periodic material (a perfect crystal), the set of such possible configurations, or the hull [START_REF] Bellissard | Coherent and dissipative transport in aperiodic solids[END_REF], is simply the periodic unit cell. Indeed, choosing the origin at points which di er only by a lattice vector results in identical configurations. This invariance by lattice translations is what enables the classical use of the Bloch theorem to reduce the Hamiltonian operator set on the whole space to a family of easily analyzed operators on the periodic unit cell, indexed by the quasimomenta k belonging to the Brillouin zone of the crystal. This is not the case for the generically incommensurate layered structures presented in Section 2.1. The change of coordinate origin (our viewpoint) is naturally associated with the action of the group R d on R d`1 by translations T that are parallel to the layers:

For a P R d , T a : # R d`1 Ñ R d`1 , x fi Ñ px 1 `a1 , . . . , x d `ad , x d`1 q . ( 2
)
Remark 2.2. Note that there is no translational symmetry in the perpendicular pd `1q stdimension, and therefore we need not include this direction in the continuous translation group. The situation is di erent in discrete (tight-binding) models, where we model hopping from layer to layer as will be seen later on.

Now we proceed with the formal definition of the hull. The positions of all the atoms are encoded in the nuclear charge distribution, a Radon measure in MpR d`1 q [2], on which the translation group R d acts naturally [START_REF] Given A P R D | space of continuous functions with compact support C c pR d`1 q through T a f pxq " f pT ´axq[END_REF]. Then the hull is the dynamical system p , R d , Tq, where is the closure of the orbit of the nuclear charge distribution measure on R d`1 generated by the atoms of all p layers under the action of R d through T.

Parameterizing this orbit is similar to describing the position of all atoms relative to the origin, given an arbitrary translation of the system. An example of this process is presented by Figure 2 (see also Figure 3 for a simple one-dimensional picture). While the initial view of all atomic positions in Figure 2a might appear quite chaotic, it can actually be e ciently encoded. Since each individual layer L j , j " 1, . . . , p is periodic, the set of all possible configurations for layer j is in one-to-one correspondence with j . Indeed, an element " j P j is the projection of the j-th layer periodic lattice on the horizontal plane, i.e., an equivalence class modulo R j , as seen in Figure 2b. The overall configuration can

1 2 3 4 h 1 = 0 h 2 h 3 h 4 b 1 b 2 b 3 b 4 1 + R 1
Figure 3: Schematic representation of a configuration Ê " p" 1 , " 2 , " 3 , " 4 q P for d " 1 and p " 4. The vector " j on the figure corresponds to the unique representative of the equivalence class " j in the periodic unit cell p j and " 1 `R1 is the set of the lattice points of the first layer in the configuration Ê. thus be parameterized as an element of

" 1 ˆ¨¨¨ˆ p . (3)
The hull therefore has the topology of a pdpq-dimensional torus. For any given configuration Ê " p" 1 , . . . , " p q, the nuclear charge distribution is

fl nuc Ê " p ÿ j"1 ¨ÿ p j P" j `Rj T p j m j p¨´h j e d`1 q '.
The support of the measure fl nuc Ê is the discrete set:

$ ' & ' % L Ê :" p § j"1 L Ê j Ä R d`1
(set of all atomic positions), with L Ê j :" Supppm j q `"j `Rj `hj e d`1

(atomic positions in layer j).

(

) 4
There exists a natural action of the additive R d group on , corresponding to an horizontal translation of the whole system in the previous parameterization, defined by

T a p" 1 , . . . , " p q " `"1 `a, . . . , " p `a˘. (5)
For convenience, we use the same notation T to denote the action of the translation group R d on R d`1 , formula (2), and on , formula [START_REF] Cao | Superlattice-induced insulating states and valley-protected orbits in twisted bilayer graphene[END_REF]. Finally, the dynamical system p , R d , Tq is equipped with a probability measure P which encodes the relative occurrences of the various configurations. This measure should be invariant by the translation action to reflect the spatial homogeneity of the system. There is here a unique such measure under the incommensurability condition given in the following definition. Let us denote the dual (or reciprocal) lattice of any cocompact [START_REF]k P R d | k ¨n P Z, @n P R ( . When R is generated by the columns of a matrix[END_REF] lattice R by R ˚.

Definition 2.3. The collection of cocompact lattices

R 1 , ¨¨¨, R p of R d is called incom- mensurate if we have for any p-tuple pk 1 , . . . , k p q P R 1 ˆ¨¨¨ˆR p , p ÿ j"1 k j " 0 ô k j " 0 @j " 1, . . . , p. (6)
This definition, surprising at first, defines incommensurability as the absence of constructive interferences, or Bragg reflections, between the lattices. The following result proves that it is also the right condition for the layered system to have a homogeneous character, i.e., all possible configurations will be visited uniformly as we translate our viewpoint along the horizontal place.

Proposition 2.4. Let R 1 , ¨¨¨, R p of R d be cocompact lattices of R d , j " R d {R j ,
lim rÑ8 1 |B r | ˆBr f pT ´aÊqda " ˆ f dP, (7)
where B r is the ball of radius r centered at zero. Proof. The first point follows from the definition (5) and the translation invariance of the Lebesgue measure. To show that incommensurability implies ergodicity, we will make use of the following lemma: Lemma 2.5. Let R be a lattice in R d and k P R d . Then:

lim rÑ8 1 # pR X B r q ÿ nPRXBr e 2ifik¨n " 0 ô k R R ˚. (8)
It is clear that if k P R ˚, then the weighted sum in the left-hand side in ( 8) is always equal to 1 for all r °0, and thus the limit is not zero. Let pc 1 , . . . , c d q be a basis of R and pc 1 , . . . , c d q the associated dual basis of R ˚(c j ¨ck " " jk ). We assume now that k R R ˚,

and expanding k " k 1 c 1 `¨¨¨`k d c d , we suppose that k 1 R Z without loss of generality. Let E K
r be the projection of R X B r onto the last pd ´1q lattice coordinates:

E K r " # pn 2 , . . . , n d q P Z d´1 ˇˇˇD n 1 P Z, d ÿ j"1 n j c j P R X B r + .
There exists for any pd ´1q-tuple r n " pn 2 , . . . , n d q in E K r two integers N ȓ pr nq such that for any n

1 P Z, ∞ d j"1 n j c j belongs to R X B r if and only if N ŕ pr nq § n 1 † N r pr nq.
We decompose accordingly the sum in (8):

ÿ nPRXBr e 2ifik¨n " ÿ r nPE K r N r pr nq´1 ÿ n 1 "N ŕ pr nq e 2ifipk 1 n 1 `r k¨r nq " ÿ r nPE K r e 2ifi r k¨r n ˜e2ifik 1 N ŕ pr nq ´e2ifik 1 N r pr nq 1 ´e2ifik 1 ¸.
We can bound the number of elements in pR X B r q from below by c 1 r d and of E K r by c 2 r d´1 where c 1 , c 2 °0 are two geometrical constants depending only on R. Thus ˇˇˇˇ1

# pR X B r q ÿ nPRXBr e 2ifik¨n ˇˇˇˇ § ˆ2 |1 ´e2ifik 1 | ˙c2 c 1 r . Since k 1 R Z, Lemma 2.
5 is proved. Birkho property. Let us now suppose that the collection of lattices R 1 , ¨¨¨, R p is incommensurate. Let f P Cp q and Ê " p" 1 , . . . , " p q P . Let Á °0 and T Á be a trigonometric polynomial such that }f ´TÁ } 8 § Á. T Á is a finite linear combination of Fourier factors of the form

G k 1 ,...,kp : p" 1 , . . . , " p q fi Ñ e 2ifipk 1 ¨"1 `¨¨¨`kp¨" p q with pk 1 , . . . , k p q P R 1 ˆ¨¨¨ˆR p . If k j " 0 for all j " 1, . . . , p, then clearly lim rÑ8 1 |B r | ˆBr G 0,...,0 pT ´aÊqda " 1 " ˆ G 0,...,0 dP.
Otherwise, we assume without loss of generality that k 2 ‰ 0, and we approximate the ball B r for r °0 by a union of translated unit cells,

B 1 r " § n 1 PR 1 XBr ´p 1 `n1 ¯where p 1 " E 1 r´1{2, 1{2q d .
Averaging over this approximate ball, we find:

1 |B 1 r | ˆB1 r G k 1 ,...,kp pT ´aÊqda " | p 1 | ´1 # pR 1 X B r q ÿ n 1 PR 1 XBr ˆp 1 e 2ifi ∞ p j"1 k j ¨p" j ´a´n 1 q da " ˆ|p 1 | ´1 ˆp 1 e 2ifi ∞ p j"1 k j ¨p" j ´aq da ˙1 # pR 1 X B r q ÿ n 1 PR 1 XBr e ´2ifi ∞ p j"2 k j ¨n1 .
We assume that the lattices R 1 , ¨¨¨, R p are incommensurate, and thus ∞ p j"2 k j P R 1 if and only if k j " 0 for all j " 2, . . . , p. This is not possible since k 2 ‰ 0, and we therefore deduce from Lemma 2.5 that, uniformly in Ê,

lim rÑ8 1 |B 1 r | ˆB1 r G k 1 ,...,kp pT ´aÊqda " 0.
Note that for some C °0 independent of r, we have |B r B 1 r | § Cr d´1 (where A B denotes the symmetric di erence of the sets A and B). Thus we conclude @pk 1 , . . . , k p q P pR 1 ˆ¨¨¨ˆR p qz t0u , lim

rÑ8 1 |B r | ˆBr G k 1 ,...,kp pT ´aÊqda " 0, uniformly in Ê.
As a result, the trigonometric polynomial T Á satisfies for r °0 large enough, uniformly in Ê, ˇˇˇ1

|B r | ˆBr T Á pT ´aÊqda ´ˆ T Á dP ˇˇˇ § Á.
Finally, since }f ´TÁ } 8 § Á we obtain the Birkho ergodic formula: for r large enough, uniformly in Ê, ˇˇˇ1

|B r | ˆBr f pT ´aÊqda ´ˆ f dP ˇˇˇ § 3Á. (9)
Ergodicity. Let B be a measurable subset of invariant under the action of R d . Let Á °0 and f Á be a continuous function such that }‰ B ´fÁ } L 1 p q § Á. By invariance of B under translations, we have for all a P R d ,

}f

Á ˝Ta ´fÁ } L 1 p q § 2Á. We can bound the di erence between f Á and its Birkho means:

› › › › 1 |B r | ˆBr f Á ˝T´a da ´fÁ › › › › L 1 p q " ˆ ˇˇˇ1 |B r | ˆBr f Á pT ´aÊqda ´fÁ pÊq ˇˇˇd PpÊq § 1 |B r | ˆBr ˆ |f Á pT ´aÊq ´fÁ pÊq| dPpÊqda § 2Á,
where we have used the triangle inequality for integrals and Fubini's theorem for nonnegative functions. Since the Birkho means of f Á converge uniformly in Ê for r Ñ 8 to

´ f Á dP, we deduce that › › › › ˆ f Á dP ´fÁ › › › › L 1 p q § 2Á.
As a consequence, we have by the triangle inequality,

}‰ B ´PpBq} L 1 p q § }‰ B ´fÁ } L 1 p q `› › › › f Á ´ˆ f Á dP › › › › L 1 p q `› › › › ˆ f Á dP ´PpBq › › › › L 1 p q § 4Á.
Since this holds for any Á °0, we conclude that ‰ B is constant a.s., and therefore PpBq P t0, 1u. Thus P is ergodic.

Non-averaging case. If the lattices R 1 , . . . , R p are not incommensurate, there exists a particular nonzero combination pk 1 , . . . , k p q P R 1 ˆ¨¨¨ˆR p such that ∞ p j"1 k j " 0. Then, we have for the corresponding Fourier factor,

G k 1 ,...,kp pT ´aÊq " G k 1 ,...,kp pÊq, @pa, Êq P R d ˆ .
The function G k 1 ,...,kp is then invariant by the R d -action and not constant. Therefore the dynamical system p , R d , T, Pq is not ergodic in this case.

Remark 2.6. When the lattices R

1 , ¨¨¨, R p are commensurate, i.e., ì p j"1 R j forms a cocompact superlattice of R d , then the orbits through the action T of R d on are no longer dense in , rather they form lower dimensional submanifolds of which are then distinct hulls for the possible configurations. In this case, Bloch theory allows one to study each of these nonequivalent configuration sets individually. Remark 2.7. Note that the lattices R 1 , ¨¨¨, R p can also be neither commensurate nor incommensurate. We do not know if one can always find a sensible parameterization of the hull of a particular configuration in this case, since its orbit under R d may not be dense in , nor has a simple geometry.

Disordered multilayer systems

The case of disordered multilayer systems is more involved. For the sake of both brevity and clarity, we consider the specific, but representative, example of a bilayer system whose bottom layer lays on a periodic substrate modeled by an external R 0 -periodic potential, and whose top layer may have defects. We denote by R 0 the periodic lattice of the substrate, and by R 1 and R 2 the periodic lattices of the bottom and top layers respectively, and we assume that R 0 , R 1 and R 2 are incommensurate in the sense of Definition 2.3. We then assume that the defects are such that -the overall geometry of the system is not modified: the periodic lattice R 2 is still appropriate to describe the configurations of the system; the di erence with the case of "perfect" homogeneous systems dealt with in the previous section, is that the nuclear distribution inside the unit cells of the top layer is not the same in each cell;

-in each cell in the top layer, the motif can be one of the following: pm Without loss of generality, we can assume that the center of the bottom layer is contained in the plane x d`1 " 0, and that the center of the top layer is contained in the plane x d`1 " h °0. For each point in r P R d , we introduce the decomposition of r associated with the lattice R 2 defined as r " rrs 2 `tru

2 with rrs 2 P R 2 and tru 2 P p 2 " E 2 r´1{2, 1{2q d . ( 10
)
For this example, the hull is the dynamical system p , R d , •, Pq where

-" 0 ˆ 1 ˆ 2 ˆt0, 1, ¨¨¨, Du R 2 , the configuration Ê " p" 0 , " 1 , " 2 , tt n 2 u n 2 PR 2 q P corresponding to the nuclear distribution fl nuc Ê " ÿ p 1 P" 1 `R1 T p 1 m 1 `ÿ p 2 P" 2 `R2 T p 2 m pt rp 2 s 2 q 2 p¨´he d`1 q.
Here t rp 2 s 2 P t0, 1, ¨¨¨, Du is the type of the motif carried by the cell of the top layer centered at p 2 `he d`1 . Note that " 0 does not appear in the definition of fl nuc Ê , since we do not take into account the geometry relaxation of the bilayer system due to the interaction with the substrate. On the other hand, the substrate generates a R 0 -periodic potential which modifies the electronic structure of the bilayer system, and this potential depends on " 0 , the relative position of the substrate with respect to the bilayer system; -• is the action of the group R d on defined by

@a P R d , @Ê " p" 0 , " 1 , " 2 , tt n 2 u n 2 PR 2 q, • a pÊq " p" 0 `a, " 1 `a, " 2 `a, t n 2 ´rt" 2 u 2 `as 2
( n 2 PR 2 q; -denoting by µ the probability on the set t0, 1, ¨¨¨, Du with law pp 0 , p 1 , ¨¨¨, p D q, the ergodic probability P on is defined as

P " P 0 b P 1 b P 2 b P dis
, where P j is the uniform probability distribution on j , and P dis " µ bR 2 is the probability distribution on the disorder.

C ˚-algebra formalism for tight-binding models

In this section, we adapt to the tight-binding modeling of perfect multilayer atomic heterostructures the C ˚-algebra formalism which was extensively used by Belissard and collaborators [START_REF] Bellissard | Coherent and dissipative transport in aperiodic solids[END_REF][START_REF] Bellissard | The noncommutative geometry of the quantum Hall e ect[END_REF] to model and analyze transport in aperiodic solids.

In the framework of tight-binding models, it is appropriate (see Remark 2.2) to use the following alternative definition of the hull:

-D " t1, ¨¨¨, pu ˆ ; -t is the action of the group G D " Z{pZ ˆRd defined by: for all a " p-, aq P G D , and all pj, " 1 , ¨¨¨, " p q P D , t a pj, "

1 , ¨¨¨, " p q " pj ´-, " 1 `a, ¨¨¨, " p `aq, where the ´sign in j ´refers to the natural action of Z{pZ on t1, ¨¨¨, pu (if -" k `pZ, j ´-" ppj ´k ´1q mod pq `1); -P D is as usual the uniform probability distribution on D , and it is ergodic if and only if the lattices R 1 , . . . , R p are incommensurate. More explicitly, for non-disordered multilayers, D consists of p copies of 1 ˆ¨¨¨ˆ p indexed by j P t1, ¨¨¨, pu such that the reference lattice site of the configuration is in layer j.

In the tight-binding representation, the quantum states are expanded on a finite number of orbitals in each periodic unit cell, for example a set of maximally localized Wannier orbitals [START_REF] Kohn | Analytic properties of Bloch waves and Wannier functions[END_REF][START_REF] Marzari | Maximally localized generalized Wannier functions for composite energy bands[END_REF][START_REF] Wannier | The structure of electronic excitation levels in insulating crystals[END_REF]. We therefore introduce, for each layer k, a finite set k of tight-binding orbitals per unit cell. For a given configuration Ê " pj, "

1 , ¨¨¨, " p q, the tight-binding orbital of type n P k associated with the lattice site of layer k located at point x `hk e d`1 P " k `Rk `hk e d`1 is indexed by the triplet ppk ´jq mod p, x, nq P Z{pZ ˆp" k `Rk q ˆ k . The integer ppk ´jq mod pq, considered as a element of Z{pZ, accounts for the vertical jump (in terms of number of layers and modulo p) to go from layer j, which contains the origin in the configuration Ê, to layer k (see also the graphical explanation in the one-dimensional case in Figure 4). This way of labeling the orbitals turns out to be well suited to the C ˚-algebra formalism introduced in the following two subsections.

The infinite set Ê of all tight-biding orbitals of the p layers is then indexed by:

Ê " p § k"1 tpk ´jq mod pu ˆp" k `Rk q ˆ k . ( 11
)
In the special case when there is only one orbital per unit cell, each k only contains one element, and can therefore be omitted in the above definition of Ê .

In the configuration Ê, the space of quantum states for the tight-binding model is

H Ê " ¸2p Ê q, (12)
and observables such as the Hamiltonian are described as linear operators on H Ê .

Abstract setting

Let us first briefly recall in this section the general formalism of groupoid C ˚-algebras.

We refer e.g. to [START_REF] Bellissard | The noncommutative geometry of the quantum Hall e ect[END_REF] for a more in-depth presentation of these mathematical objects. We further simplify the presentation by assuming that there is only one tight-binding orbital per unit cell in each layer, so that

Ê " p § k"1 tpk ´jq mod pu ˆp" k `Rk q Ä G D . ( 13
)
Note that this is not a restriction of the formalism (see Remark 3.1 below). The first step is to construct a groupoid based on the canonical transversal X of the hull D : X " Ê P D | Ê " pj, " 1 , . . . , " p q; " j " 0

( . ( 14
)
The idea is that while D indexes the possible viewpoints from any point in the layer planes, each element of X represents a possible unique viewpoint from the position of a lattice site, which is then chosen as the origin. This is a more appropriate approach in the case of tight-binding models. Associated with the transversal X is the groupoid pXq defined as follows:

pXq " tpÊ, aq P X ˆGD | t ´aÊ P Xu . [START_REF] Marzari | Maximally localized generalized Wannier functions for composite energy bands[END_REF] For Ê " pj, " 1 , ¨¨¨, " p q P D and a " p-, aq P G D , pÊ, aq P pXq if and only if " j " 0 and " j`-" a. Geometrically, pXq indexes all the possible jumps between two lattice sites in all possible configurations of the multilayer system. The physical space vector of the jump associated with pÊ, aq is A " a `ph j`-´hj qe d`1 . The groupoid is equipped with a set of three operations, respectively the range r : pXq Ñ X, source s : pXq Ñ X, and composition ˝: pXq ˆ pXq Ñ pXq, satisfying: rpÊ, aq " Ê, spÊ, aq " t ´aÊ, pÊ, aq ˝pt ´aÊ, bq " pÊ, a `bq.

The fiber pÊq is defined as r ´1ptÊuq for any Ê P X.

In a second step, we define the ˚-algebra A 0 of continuous functions with compact support defined on pXq with values in C, endowed with the following composition laws and ˚-operator: p⁄f qpÊ, aq " ⁄f pÊ, aq, pf `gqpÊ, aq " f pÊ, aq `gpÊ, aq, pf ˚gqpÊ, aq " ÿ pÊ,bqP pÊq f pÊ, bqgpt ´bÊ, a ´bq,

f ˚pÊ, aq " f pt ´aÊ, ´aq. (17)
The ˚-algebra A 0 has an identity, denoted by 1: 1pÊ, aq " " a,0 .

Remark 3.1. In the case when there are several atoms per unit cell, and/or each atom carries more than one orbital, the functions f in [START_REF] Given A P R D | space of continuous functions with compact support C c pR d`1 q through T a f pxq " f pT ´axq[END_REF] do not take their values in C but f pÊ, aq is in C N j ˆNj`-when Ê " pj, " 1 , . . . , " p q and a " p-, aq, where N k " #p k q is the total number of tight-binding orbitals carried by the atoms in the unit cell of layer k. Products of the form f pÊ, bqgpt ´bÊ, a ´bq should then be understood as matrix products.

The ˚-algebra A 0 can be mapped onto a ˚-algebra of bounded linear operators acting on the space of quantum states H Ê " ¸2p Ê q, see ( 12), via the representation formula: fi Ê pf q"pxq " ÿ pÊ,yqP pÊq f pt ´xÊ, y ´xq "pyq, @" P H Ê , @x P Ê . (

In the above formula, we have implicitly used the fact that the set Ê defined by ( 13) coincides with y P G D | pÊ, yq P pÊq ( .

Configuration ! 0 Y X x y x = ( 1, x) y = (1, y) 0 Y X y x y x = (2, y x) Configuration t x !
Figure 4: Graphical explanation of the representation formula (18) for d " 1, p " 4, and for a configuration Ê " p3, " 1 , " 2 , 0, " 4 q P D . The points x " p2 ´3, xq " p´1, xq and y " p4 ´3, yq " p1, yq of Ê are combined to form an element y ´x " p2, y ´xq such that pt ´xÊ, y ´xq P pt ´xÊq .

The following covariance condition holds: for a given " " pÊ, aq P pXq, fi t ´aÊ pf q " Tp"q ´1fi Ê pf qTp"q,

where the translations Tp"q : H t ´aÊ Ñ H Ê are unitary operators defined by Tp"q"pxq " "pt ´axq, @" P H t ´aÊ , @x P Ê . (

) 20
Remark 3.2. When d " 2, magnetic fields can be included in the description through a Peierls substitution term as follows. We assume that there is only one orbital per unit cell, and that the orbitals are carried by atoms located at the lattice sites. Let B be an antisymmetric 3 ˆ3 matrix representing the magnetic field. For x " pl, xq and y " pm, yq in Ê with Ê " pj, " 1 , . . . , " p q, we define the corresponding Peierls phase term as

Ê B px, yq " fi 3 ÿ ‹,µ"1 B ‹,µ X ‹ Y µ , ( 21
)
where X " x `ph j`l ´hj qe 3 and Y " y `ph j`m ´hj qe 3 are the physical space vectors of the jumps from the origin to the sites x and y.

The product and representation formulae are then modified as pf ˚B gqpÊ, aq " ÿ pÊ,bqP pÊq f pÊ, bqgpt ´bÊ, a ´bqe i Ê B pa,bq , fi Ê pf q"pxq " ÿ pÊ,yqP pÊq f pt ´xÊ, y ´xq e ´i Ê B px,yq "pyq, @" P H

Ê , x P Ê , ( 22
)
where we again used the fact that Ê " y P G D | pÊ, yq P pÊq ( . In this case, the following covariance condition holds: for a given " " pÊ, aq P pXq, fi t ´aÊ pf q " U p"q ´1fi Ê pf qU p"q, (23

)
where the magnetic translations U p"q : H t ´aÊ Ñ H Ê are unitary operators defined by U p"q"pxq " exp ˆi ˆrX´A,Xs A pyq ¨dy ˙"pt ´axq, @" P H t ´aÊ , @x P Ê . (

) 24
Here, A is a vector potential giving rise to the magnetic field B, and rX ´A, Xs is the line segment joining X ´A to X in R 3 .

Note that fi

Ê pf q is hermitian if f " f ˚. This representation induces a C ˚norm on A 0 , defined by: }f } " sup ÊPX }fi Ê pf q}, where the norm on the right hand side is the operator norm on LpH Ê q, the space of bounded linear operators on H Ê . We then construct in a third and final step the C ˚algebra A as the completion of A 0 for this norm.

Integro-di erential calculus

Integration and derivations can be constructed on the C ˚-algebra as follows. Let P be a probability measure on X invariant by pXq-action. A formal integration is then obtained on A 0 as follows: if f P A 0 ,

T P pf q " ˆX dPpÊqf pÊ, 0q. ( 25
)
This trace is positive, T P pf ˚˚f q • 0. It is faithful (T P pf ˚˚f q " 0 i f " 0) when the support of the measure P is X.

Remark 3.3. This trace coincides with the trace per unit volume in R d of the corresponding operator fi Ê pf q thanks to the Birkho property [START_REF] Cazeaux | Analysis of rippling in incommensurate one-dimensional coupled chains[END_REF].

One then defines L p -norms by setting for 1

§ p † 8, }f } p " T P `pf ˚f ˚qp{2 ˘1{p @f P A 0 . ( 26
)
Closure of A 0 with respect to the L p -norms defines Banach spaces L p pA, T P q. Derivations are introduced next by setting

B j f pÊ, xq " ix j f pÊ, xq, f P A 0 . ( 27
)
Given a multi-index -" p-1 , . . . ,d q we will use the notation

B -f for B -1 1 . . . B - d d .
The derivations defined above satisfy the fundamental properties of derivation operators: they commute, they are ˚-derivations in the sense that B j pf ˚q " pB j f q ˚, and satisfy the Leibniz rule B j pf ˚gq " B j f ˚g `f ˚Bj g . Finally one has the operator representation:

fi Ê pB j f q " ´i " x Ê j , fi Ê f ‰ , ( 28
)
where x Ê " px Ê 1 , . . . , x Ê d q is the position operator on L Ê .

Remark 3.4. In the periodic case, these derivations correspond (by Fourier transform) to derivation in quasi-momentum space.

Analytic spectral calculus A basic notion in operator algebras is that of the resolvent sets, the spectrum, and the introduction of a spectral calculus by a complex contour integral. Since the operators we are considering here are bounded, this is straightforward and we refer e.g. to [START_REF] Arveson | A short course on spectral theory[END_REF] for the details. For a given element f P A, its resolvent set flpf q and spectrum ‡pf q are flpf q " tz P C s.t. pz1 ´f q is invertibleu, ‡pf q " Czflpf q.

Furthermore, the resolvent set is open in C, and the resolvent function z fi Ñ pz1 ´f q ´1 is an algebra-valued analytic function of z P flpf q. An analytical calculus can then be defined on A as follows. Let F ‡pf q by the algebra of C-valued functions which are analytic in an open neighborhood of ‡pf q. Then we define a homomorphism of algebras

F ‡pf q Q " fi Ñ 1 2ifi ˛C "pzqpz1 ´f q ´1dz P A, ( 29
)
where C is a contour surrounding ‡pf q in the analyticity domain of ". The right-hand side of this mapping is independent of the particular choice of contour, and usually noted as "pf q.

Perfect bilayers

Let us identify further the groupoid and the associated C ˚-algebra in the specific geometry of a bilayer system with no disorder.

Transversal and groupoid. We start by identifying the transversal in the case of the hull D of a perfect incommensurate bilayer. We have

D " t1, 2u ˆ 1 ˆ 2 ,
and X " tp1, 0, " 2 q ; " 2 P 2 u Y tp2, " 1 , 0q ; " 1 P 1 u . This leads to the identification:

X " X 1 Y X 2 where # X 1 " 2 , X 2 " 1 . ( 30
)
The set X 1 (resp. X 2 ) describes the possible configurations from the point of view of a given lattice site of layer L 1 (resp. L 2 ). In a configuration " 2 P X 1 , which corresponds to Ê " p1, 0, " 2 q, the lattice sites of L 1 are located at R 1 while the lattice sites of L 2 are located at " 2 `R2 `he d`1 , where h is the distance between the two layers. In a configuration " 1 P X 2 , which corresponds to Ê " p2, " 1 , 0q, the lattice sites of L 2 are located at R 2 while the lattice sites of L 1 are located at " 1 `R1 ´he d`1 . The decomposition (30) of the transversal allows us to identify a block decomposition of the groupoid:

pXq " ˛ 11 Y ˛ 12 Y ˛ 21 Y ˛ 22 " ˛ pX 1 , X 2 q, ( 31
)
where the set of arrows ˛ 11 and ˛ 22 include all intralayer jumps:

# ˛ 11 " tp" 2 , mq ; " 2 P X 1 , m P R 1 u , ˛ 22 " tp" 1 , nq ; " 1 P X 2 , n P R 2 u ,
while ˛ 21 describes jumps from the first to the second layer, and ˛ 12 jumps from the second to the first layer. Let us note that the above representation would be redundant in the case of interlayer jumps, e.g., for p" 2 , qq in ˛ 12 . Indeed, " 2 P X 1 can be deduced from q P " 2 `R2 as its equivalence class modulo R 2 . We will thus denote these arrows solely as arrow vectors q, observing that these can take their values in all of R d :

# ˛ 12 " q; q P R d ( , ˛ 21 " p; p P R d ( . ( 32
)
This allows us to make explicit the three groupoid operations, namely the range map

r : ˛ Ñ ˛ , $ ' ' ' ' & ' ' ' ' % ˛ 11 Ñ X 1 , p" 2 , mq fi Ñ " 2 , ˛ 12 Ñ X 1 , q fi Ñ q R 2 , ˛ 21 Ñ X 2 , p fi Ñ p R 1 , ˛ 22 Ñ X 2 , p" 1 , nq fi Ñ " 1 , , (33)
the source map:

s : ˛ Ñ ˛ , $ ' ' ' ' & ' ' ' ' % ˛ 11 Ñ X 1 , p" 2 , mq fi Ñ " 2 ´mR 2 , ˛ 12 Ñ X 2 , q fi Ñ ´qR 1 , ˛ 21 Ñ X 1 , p fi Ñ ´pR 2 , ˛ 22 Ñ X 2 , p" 1 , nq fi Ñ " 1 ´nR 1 , (34)
and the composition map: $ ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' % For p" 2 , mq P ˛ 11 , m 1 P R 1 : p" 2 , mq ˝p" 2 ´mR 2 , m 1 q " p" 2 , m `m1 q P ˛ 11 , for q P ˛ 12 , p `q P R 1 : q ˝p " pq R 2 , q `pq P ˛

q ˝p´q R 1 , nq " ›››Ñ q `n P ˛ 12 , ( 35
)
where we have introduced the notation r R j " r`R j , j " 1 . . . p to denote the equivalence classes in j of a given point r P R d . Finally, the fiber pÊq can be determined: ˛ p" 2 q " p1,0," 2 q " L p1,0," 2 q , for " 2 P X 1 , (36a) ˛ p" 1 q " p2," 1 ,0q " L p2," 1 ,0q , for

" 1 P X 2 . ( 36b
)
where L Ê is the set of points defined by (4). We obtain for the groupoid the discrete counterpart to the ergodic property of the continuous hull, Proposition 2.4. Let d" denote the usual Lebesgue measure on X

Proposition 3.5. Let P be the probability measure on X with uniform density p| 1 | | 2 |q ´1d". 1. P is invariant by the groupoid action.

2. The dynamical system pX, pXq, t, Pq is uniquely ergodic if and only if the lattices R 1 , ¨¨¨, R p are incommensurate. In this case, we have the Birkho property: for any f P CpXq and Ê P X,

lim rÑ8 1 # pB r X L Ê q ÿ aPBrXL Ê f pt ´aÊq " ˆX f dP, ( 37
)
where B r is the ball of radius r centered at the origin. The proof of this result follows the same lines as the proof of Proposition 2.4. We therefore omit it for the sake of brevity.

Bilayer C ˚-algebra. For the case of periodic incommensurate bilayers with one orbital per unit cell, we can give a comprehensive description of the abstract algebra defined above. First, elements of H Ê " ¸2p Ê q can be seen as: -elements of H " 2 " ¸2pR 1 q ' ¸2p" 2 `R2 q for Ê " " 2 P X 1 , -or elements of H " 1 " ¸2p"

1 `R1 q ' ¸2pR 2 q for Ê " " 1 P X 2 .

Then, given the decomposition (31) of pXq, it makes sense to write a block decomposition of functions f P C ˚p pXqq as:

f " " f 11 f 12 f 21 f 22 ⇢ (38)
where

f 11 : ˛ 11 " 2 ˆR1 Ñ C, f 12 : ˛ 12 " R d Ñ C, f 21 : ˛ 21 " R d Ñ C, f 22 : ˛ 22 " 1 ˆR2 Ñ C. (39)
Note that the decomposition (38) of the tight-binding hopping parameters into intra-and inter-layer terms is naturally the parameterization used in physics [START_REF] Fang | Ab initio tight-binding hamiltonian for transition metal dichalcogenides[END_REF]. In particular, inter-layer coe cients are usually represented directly as a continuous function of the relative position of the atoms as in (39), see e.g., [START_REF] Fang | Electronic structure theory of weakly interacting bilayers[END_REF]. Let us now write the ˚product defining the algebra:

-for " 2 P X 1 and m P R 1 , pf ˚gq 11 p" 2 , mq "

ÿ m 1 PR 1 f 11 p" 2 , m 1 qg 11 p" 2 ´m1 R 2 , m ´m1 q `ÿ q 1 P" 2 `R2 f 12 pq 1 qg 21 pm ´q1 q; (40a)
-for q P R 2 , pf ˚gq 12 pqq "

ÿ 1 PR 1 f 11 pq R 2 , m 1 qg 12 pq ´m1 q `ÿ n 1 PR 2 f 12 pq ´n1 qg 22 pn 1 ´qR 1 , n 1 q; (40b)
-for p P R 2 , pf ˚gq 21 ppq "

ÿ n 1 PR 2 f 22 pp R 1 , n 1 qg 21 pp ´n1 q `ÿ m 1 PR 1 f 21 pp ´m1 qg 11 pm 1 ´pR 2 , ´m1 q; (40c) -for " 1 P X 2 and n P R 2 , pf ˚gq 22 p" 1 , nq " ÿ n 1 PR 2 f 22 p" 1 , n 1 qg 22 p" 1 ´n1 R 1 , n ´n1 q `ÿ p 1 P" 1 `R1 f 21 pp 1 qg 12 pn ´p1 q. (40d)
We can also write the ˚operation as follows: for " 2 P X

Finally, the representation formula writes as follows:

1. for " 2 P X 1 and " " p" 1 , " 2 q P ¸2pR 1 q ' ¸2p" 2 `R2 q, fi " 2 pf q" can be decomposed in ¸2pR

1 q ' ¸2p" 2 `R2 q as $ ' ' ' ' ' & ' ' ' ' ' % `fi" 2 pf q" ˘1 pmq " ÿ m 1 PR 1 f 11 `"2 ´mR 2 , m 1 ´m˘" pm 1 q `ÿ q 1 P" 2 `R2 f 12 pq 1 ´mq "pq 1 q, `fi" 2 pf q" ˘2 pqq " ÿ

m 1 PR 1 f 21 pm 1 ´qq "pm 1 q `ÿ n 1 PR 2 f 22 `´q R 1 , n 1 ˘"pq `n1 q,
for all m P R 1 and q P " 2 `R2 ; (42) 2. for " 1 P X 2 and " " p" 1 , " 2 q P ¸2p" 1 `R1 q ' ¸2pR 2 q, fi " 1 pf q" can be decomposed in ¸2p"

1 `R1 q ' ¸2pR 2 q as $ ' ' ' ' ' & ' ' ' ' ' % `fi" 1 pf q" ˘1 ppq " ÿ

m 1 PR 1 f 11 `´p R 2 , m 1 ˘"pp `m1 q `ÿ n 1 PR 2 f
12 pn 1 ´pq "pn 1 q, `fi" 1 pf q" ˘2 pnq " ÿ

p 1 P" 1 `R1 f 21 pp 1 ´nq "pp 1 q `ÿ n 1 PR 2 f 22 `"1 ´nR 1 , n 1 ´n˘" pn 1 q,
for all p P " 1 `R1 and n P R 2 .

(43)

Trace per unit volume and integro-di erential calculus. A classical consequence [START_REF] Bellissard | The noncommutative geometry of the quantum Hall e ect[END_REF] of the ergodicity of P on X under the action of pXq, Proposition 3.5, is that we can characterize uniquely P under the condition that it is a trace per unit volume: Proposition 3.6. When R 1 and R 2 are incommensurate, the invariant, ergodic probability measure P is uniquely defined as a trace per unit volume in the sense that it satisfies: for any f P C ˚p pXqq and Ê P X,

T P pf q " lim rÑ8 1 # pB r X L Ê q Tr pfi Ê pf q| Br q . ( 44

Noncommutative Kubo formula

Electronic transport can be modeled as the result of the interplay between the quantum evolution in the presence of a uniform electric field and dissipation mechanisms such as scattering events which depend on the environment of the charge carrier (electron or hole) with charge q. As discussed in the lecture notes [START_REF] Bellissard | Coherent and dissipative transport in aperiodic solids[END_REF], models for these dissipation mechanisms can be obtained by considering microscopic, many-body electron systems coupled with an environment such as a phonon bath, and then integrating out the degrees of freedom of the environment to obtain an e ective single-particle model. As in the Drude model, dissipation mechanisms are represented by discrete scattering events (collisions) with Poisson-distributed independent time delays between successive collisions. Given any one-particle density matrix fl as an initial state (i.e. a positive element of the C ˚-algebra associated with a quantum system), the e ective, collision-averaged one-particle time evolution is given by the Liouville equation [START_REF] Schulz-Baldes | A kinetic theory for quantum transport in aperiodic media[END_REF]

dfl dt `Lh´qEptq¨x pflq " ´1 ´p Ÿ • pflq, (46)
where h is the e ective single-electron tight-binding Hamiltonian, -E is a constant or time-harmonic spatially uniform electric field with frequency p Ê, -x is the position operator and 1 is the identity operator, -q is the charge of the carrier (hole or electron), -L h´qEptq¨x " i{~rh, ¨s´q ~E ptq¨Ò is the Liouvillian operator governing the one-particle time evolution in the absence of collisions, p Ÿ ˚is a scattering-event-averaged collision e ciency operator, -• is the mean collision time of a Poisson process with law e ´t{• dt{• governing the time delay between the independent scattering events (collisions).

The linear conductivity tensor is then readily available in the framework of C ˚-algebras by the famous noncommutative Kubo formula [START_REF] Schulz-Baldes | A kinetic theory for quantum transport in aperiodic media[END_REF]: ‡ ij pp Êq "

q 2 ~2 T P `Bi h rp1 ´p Ÿ ˚q{• `Lh ´ip Ês ´1 B j f -,µ phq ˘, (47)
where f -,µ phq " 1 1`e -ph´µq is the Fermi-Dirac one-particle density matrix in the grandcanonical equilibrium with chemical potential (Fermi level) µ and at temperature T such that -" 1{k B T . In the relaxation time approximation (RTA), the e ective relaxation operator p1 ṕ Ÿ ˚q{• is replaced by 1{• rel , a single relaxation time which depends in general strongly on the temperature. Full details of the modeling assumption and derivation of (47) can be found in the references [START_REF] Bellissard | Coherent and dissipative transport in aperiodic solids[END_REF][START_REF] Bellissard | The noncommutative geometry of the quantum Hall e ect[END_REF][START_REF] Prodan | Quantum transport in disordered systems under magnetic fields: A study based on operator algebras[END_REF][START_REF] Schulz-Baldes | A kinetic theory for quantum transport in aperiodic media[END_REF].

Numerical example

We propose in this section a minimalistic one-dimensional toy model to study incommensurability e ects in multilayer systems, which is new up to our knowledge. We then present a numerical strategy and numerical results for the computation of the density of states and the conductivity. In this paper, we will restrict ourselves to scanning rational values of the lattice ratio parameter. We therefore use the traditional approach of constructing periodic supercells [START_REF] Prodan | Quantum transport in disordered systems under magnetic fields: A study based on operator algebras[END_REF], and we never actually compute directly quantities of interest in the case of incommensurate lattice ratios. In future publications, we will propose strategies based directly on the C ˚-algebra representation, possibly addressing directly the incommensurate case, and with rigorous error control. This substantial additional e ort is currently ongoing.

Description of the model

Let us consider two parallel one-dimensional crystals with lattice constants ¸1 and ¸2 normalized such that ¸1¸2 " 1.

Following the notation introduced in Sections 2 and 3, we set the lattices R j " ¸jZ and the unit cells j " R{R j " ¸jT for j " 1, 2, where T :" R{Z is the one-dimensional 1-periodic torus. The set of all possible configurations is parameterized by the hull " 1 ˆ 2 . We then form a quantum lattice model with one orbital per unit cell in each layer. The relevant parameterization is provided by the transversal X " X 1 Y X 2 with X 1 " ¸2T and X 2 " ¸1T. (49)

The corresponding Hamiltonian of the model H Ê " fi Ê phq, given by (42) and sketched on Figure 5, features two intra-chain first-neighbor hopping models with amplitude normalized to 1 for each chain, and an inter-chain coupling term with a Gaussian profile Remark 4.1. Note that knowledge of the current-current correlation measure, which depends only on the Hamiltonian, is enough to compute the conductivity for any values of the Fermi level µ or Boltzmann factor -" 1{k B T .

Kernel polynomial method

The second ingredient in our calculations is a Chebyshev polynomial expansion, commonly called the Kernel Polynomial Method [START_REF] Weiße | The kernel polynomial method[END_REF]. The idea is to rescale the Hamiltonian, H - 0 Ñ Ĥ-0 " pH - 0 ´bq{a with a, b well chosen so that the spectrum ‡p Ĥ-0 q lies in the energy range p´1, 1q.

Density of States.

The moments of the spectral measure dµ -can be computed on the basis of Chebyshev polynomials of the first kind tT m u m•0 , which form an orthonormal basis of L 2 pr´1, 1sq with respect to the weight function wpxq " 1{pfi ? 1 ´x2 q. These polynomials obey the recursion relation,

T 0 pxq " 1, T 1 pxq " x, T m`1 pxq " 2xT m pxq ´Tm´1 pxq for m • 1. (55)
From (51), we thus obtain that

µ - m " ˆR T m ˆE ´b a ˙dµ -pEq « 1 N Tr ´Tm p Ĥ- 0 q ¯" 1 N N ÿ j"1 T m p ⁄- j q, (56)
where t ⁄j u 1 §j §N is the vector of eigenvalues of Ĥ-0 . Utilizing the recursion (55), the moments pµ - m q 0 §m §M can be computed e ciently up to some fixed polynomial degree M . Now, assuming that h -has an absolutely continuous spectrum, the density of states is continuous with respect to the Lebesgue measure, dµ -pEq " ‹ -pEqdE, we can reconstruct accurately the spectral density ‹ -from the Chebyshev moments [START_REF] Weiße | The kernel polynomial method[END_REF]:

‹ -pEq " 1 fi a a 2 ´pE ´bq 2 ˜µ0 `2 M ÿ m"1 µ - m g M m T m ˆE ´b a ˙¸, (57)
where g M m " rpM ´m`1q cosp fim M `1 q`sinp fim M `1 q cotp fi M `1 qs{pM `1q are the Jackson damping coe cients designed to avoid spurious Gibbs oscillations.

Finally, we note that the values of ‹ -at the particular set of points

x k " a cos ˆfipk `1{2q M ˙`b with k " 0, . . . , M ´1, (58)
coinciding with the abscissas of the Chebyshev-Gauss numerical integration points, can be obtained through a fast cosine transform in OpM log M q operations since In conclusion, formulae (57) and (60) show that the Chebyshev expansion leads to an e cient representation of the spectral measures, relying on fast discrete cosine transforms to compute weights " - k and - kl that depend only on the matrix H - 0 . The knowledge of these coe cients enables accurate calculations of the quantities of interest for any choice of the energy E for the density of states, or of the four parameters µ, -, • rel or p Ê for the conductivity.

" - k " fi
Remark 4.2. In this work, we present some results for relatively small Hamiltonian matrices, for which a full diagonalization is feasible and provides an e ective way of computing the quantities of interest. For larger matrices, a promising alternative is the stochastic evaluation of traces which can be used to compute Chebyshev moments such as (59) needed in the Kernel Polynomial Method, losing some accuracy but enabling much larger calculations [START_REF] Weiße | The kernel polynomial method[END_REF].

Numerical results

For the purpose of the numerical examples presented in this section, we use the set of numerical parameters presented in Table 1. The numerical strategy proposed above was implemented in Julia [START_REF] Bezanson | Julia: A fast dynamic language for technical computing[END_REF]. We plot first the density of states as a function of lattice constants ratioand energy E in Figure 6. A clear fractal pattern of band gaps emerges, with continuous dependence on the lattice constants ratio parameter. Note the divergence of the density of states at the edges of the gaps due to the one-dimensional nature of the system.

The overall pattern is reminiscent of the Hofstadter butterfly [START_REF] Hofstadter | Energy levels and wave functions of bloch electrons in rational and irrational magnetic fields[END_REF], which is a paradigm of fractal structure in the density of states of an electronic Hamiltonian induced by the interplay between two length scales (lattice and magnetic field), measured by the magnetic flux through the unit cell. In particular, around -" ¸2{¸1 " 1 a large number of gaps open at the top and bottom of the spectrum. Although the resolution in this region is not very good, the similarity with Landau levels for which the energy is proportional to the magnetic field (the incommensurability parameter) is striking.

Next, we plot the conductivity as a function of lattice constants ratioand Fermi level µ in Figure 7. The same fractal pattern emerges, however a striking di erence is that the conductivity drops at the edges of the gaps (where the density of states is maximum). Strong insulating gaps occur along the spectral gaps as expected, forming a strong fan structure around -" 1.

This fan structure is repeated at a number of values ofcorresponding to rational numbers: 1{5 and 5, 1{3 and 3 are particularly strong, with weaker features at 1{4 and 4, 1{2 and 2, 3{5 and 5{3. In general, features such as gaps opening or closing appear at rational values of -. Finally, we plot in Figure 8 the conductivity again, which we plot this time as a function of the integrated density of states n -pEq " ˆE ´8 dµ -, where dµ -is the density of states measure defined by (51). The function E fi Ñ n -pEq is increasing from 0 to 1, but stays constant in a gap of the spectrum, thus in this representation the size of the insulating gaps is narrowed. Note that in experiments, the control is over the number of electrons per unit cell through gating or doping [START_REF] Sugawara | Fabrication of Li-intercalated bilayer graphene[END_REF], corresponding to the integrated density of states, and not directly over the Fermi level in general.

Surprisingly, the gaps are still clearly visible in this rescaled presentation, but appear as straight lines. This feature is reminiscent of experimental images of magneto-transport data in small twist angle bilayer graphene [START_REF] Cao | Superlattice-induced insulating states and valley-protected orbits in twisted bilayer graphene[END_REF], leading to exciting perspectives for the application of the framework presented in this paper to more realistic 2D multilayer systems.

Figure 1 :

 1 Figure 1: Example of a heterostructure: stack of five monolayers

  (a) Top view of the crystalline structure

O

  Projection of the origin into the unit cells

Figure 2 :

 2 Figure 2: Illustration of the process for labeling the local configurations of heterostructures. The vertices of the parallelograms are sites of the shifted lattices " j `Rj .

p0q 2 , m p1q 2 ,

 22 ¨¨¨, m pDq 2 q, where m p0q 2 is the periodic motif of the top layer in the absence of defects, and m p1q 2 , ¨¨¨, m pDq 2correspond to the di erent kinds of defects that can be observed. For instance, in the case of a graphene layer for which each carbon atom can adsorb a hydrogen atom, there are two carbon atoms in each cell, m p0q 2 corresponds to the case when no hydrogen atom is adsorbed, m p1q 2 (respectively, m p2q 2 ) to the case when only the first (respectively, second) carbon atom has adsorbed a hydrogen atom, and m p3q 2 to the case when the two carbon atoms have adsorbed hydrogen atoms; -the defects are independently and identically distributed in the cells of the top layer. We denote by pp 0 , p 1 , ¨¨¨, p D q (with p j °0 and

a a 2

 2 ´px k ´bq 2 ‹ -px k q "This further diminishes the numerical cost of evaluating (57).where the - kl are the raw output of the 2D fast cosine transform -

  1000 .5 .25 250 250 0

Figure 6 :

 6 Figure 6: Density of states in color scale as a function of energy and lattice constants ratios.

Figure 7 :

 7 Figure 7: Conductivity in color scale as a function of Fermi level and lattice constants ratios, with parameter choice W " .5, r 0 " .25, -" • " 250.

  and " 1 ˆ¨¨¨ˆ p , endowed with the uniform probability measure P. Then, 1. P is invariant by the translation group R d ; 2. the dynamical system p , R d , T, Pq is uniquely ergodic if and only if the lattices R

1 , ¨¨¨, R p are incommensurate. In this case, we have the Birkho property: for any f P Cp q and Ê P ,

  is bounded on A. Hence, an element f belongs to C N if and only if }B -f } † 8 for all multi-indexessuch that ∞ The extension to multilayer systems of the formalism introduced in the previous sections is straightforward, though somewhat cumbersome.

					)
	Moreover, for j " 1, ¨¨¨, d, the derivation defined by (27) extends easily: for " and n P R 1 , q P " 2 , p P " 1 , m P R 2 ,	1 2 P X
	pB j f q 11 p" pB j f q 21 ppq " ip 2 , nq " in	j f 11 p" j f 21 ppq, 2 , nq,	pB j f q 12 pqq " iq j f 12 pqq, j f q pB 22 p" 1 , mq " im j f 22 p" 1 , mq.	(45a)
	There exists also a derivation for j " d	`1
	pB d`1 f q 11 p" pB d`1 f q 21 ppq " ´if 2 , nq " 0, 21 ppq,	pB d`1 f q 12 pqq " `if 12 pqq, d`1 f q pB 22 p" 1 , mq " 0.	(45b)
	Note that the derivation operator B d`1 d j"1 -

j § N . Remark 3.7.

  Figure 5: Sketch of the geometry and tight-binding hoppings of our 1D toy model. We consider the model Hamiltonian generated by an element h of the C ˚-algebra presented in Section 3.2 and which takes the values

			2 (n 1)	2 n	2 (n + 1)
			2
			1
				1 m	1 (m + 1)
	h 11 p"	2 , mq "	# 1 if |m| " ¸1, 0 otherwise,		for "	2 P ¸2T, m P ¸1Z,
	h 12 ppq " W e ´1 2 ´p r 0 ¯2 , h 21 pqq " W e ´1 2 ´q r 0 ¯2 , h 22 p" 1 , nq " # 1 if |n| " ¸2, 0 otherwise,		for "	for p P R, for q P R, 1 P ¸1T, n P ¸2Z.

Table 1 :

 1 Choice of numerical parameters

Acknowledgments

This work was supported in part by ARO MURI Award W911NF-14-1-0247. EC and PC are grateful to Jean Bellissard for useful discussions.

depending on the distance between the lattice sites. The two parameters of the model are the maximum amplitude W and the characteristic length r 0 of the inter-chain hopping interaction terms.

Numerical approach

Using this toy model, we aim to showcase the expected e ects of continuously varying the lattice constants ratio -" ¸2{¸1 in a range around the periodic case of matched chains, ¸1 " ¸2 i.e. -" 1.

Periodic supercells

We use here the well-known numerical strategy of scanning all periodic approximations with a given total number of sites N in this range [START_REF] Hofstadter | Energy levels and wave functions of bloch electrons in rational and irrational magnetic fields[END_REF][START_REF] Prodan | Quantum transport in disordered systems under magnetic fields: A study based on operator algebras[END_REF], thus creating large periodic supercells. Allocating p atoms to the top chain and q atoms to the bottom chain so that p `q " N , we set

Note that the periodic supercell length is ? pq " p¸1 " q¸2. We then scan the range of ratios 1{6 § - § 6 by varying p from p min "

\ . The corresponding tight-binding hamiltonian matrix H - 0 with periodic boundary conditions is then assembled from the C ˚-algebra element h -(49) using a cut-o distance of 6r 0 for the Gaussian interlayer term, and choosing the configuration "

1 " " 2 " 0. Our two quantities of interest are as follows.

-the density of states, i.e., the spectral measure dµ -pEq on R with support on ‡phq is defined by ˆR "pEqdµ -pEq "

where the test function " is analytic in an open neighborhood of ‡ph -q;

-the conductivity at zero frequency given by the Kubo formula (47), which can be computed more e ciently by introducing the current-current correlation measure dM -:

The spectral measure dM -on R 2 is defined by [START_REF] Prodan | Mapping the current-current correlation function near a quantum critical point[END_REF][START_REF] Schulz-Baldes | A kinetic theory for quantum transport in aperiodic media[END_REF]:

where the test functions " 1 , " 2 are analytic in an open neighborhood of ‡ph -q. Since x fi Ñ x is not ? pq-periodic, we have used in (53) the approximate periodic di erential calculus introduced by Prodan [START_REF] Prodan | Quantum transport in disordered systems under magnetic fields: A study based on operator algebras[END_REF], such that the approximate derivation

where x fi Ñ r Xpxq is ? pq-periodic and approximates the identity x fi Ñ x near zero.

Conductivity.

Let us define the two-dimensional Chebyshev moments of the currentcurrent correlation measure dM -following (53):

To simplify the computation, let us introduce the diagonal matrix of eigenvalues D-" Diagp⁄ - 1 , . . . , ⁄ - N q, a unitary matrix of eigenvectors V -of Ĥ-0 , and the Hermitian matrix

Then,

The moments M - mn can thus be e ciently computed up to the partial degree M , using the recursion (55) as before. They can then be used to evaluate the conductivity by Chebyshev-Gauss numerical integration in (52) as follows. Denoting -,µ,• rel ,p Ê pE, E 1 q the integrand in the right-hand side of (52), we introduce approximate Chebyshev moments mn -,µ,• rel ,p Ê computed with Chebyshev-Gauss integration on the points x k defined by (58):

« p1 `"m,0 qp1 `"n,0 q ¨r´1,1s 2 -,µ,• rel ,p Ê pb `ax, b `ayqT m pxqT n pyq fi 2 ? 1 ´x2 ? 1 ´y2 dxdy, so that the function -,µ,• rel is well approximated by the expansion

Thanks to (52), the conductivity is then given by: ‡ "

Finally, we use the relations T m px k q " cos ´fimpk`1{2q M ¯and T n px l q " cos ´finpl`1{2q M ¯and we exchange the order of summations. The numerical approximation to the conductivity is then given by the quadrature formula


E P Kontar

J E Perez

L K Harra

A A Kuznetsov

A G Emslie

N L S Jeffrey

N H Bian

& B R Dennis

Turbulent kinetic energy in the energy balance of a solar flare

Keywords: numbers: 96.60.qe, 52.35.Vd, 52.65.Cc, 96.60.Iv

The energy released in solar flares derives from a reconfiguration of magnetic fields to a lower energy state, and is manifested in several forms, including bulk kinetic energy of the coronal mass ejection, acceleration of electrons and ions, and enhanced thermal energy that is ultimately radiated away across the electromagnetic spectrum from optical to X-rays. Using an unprecedented set of coordinated observations, from a suite of instruments, we here report on a hitherto largely overlooked energy component -the kinetic energy associated with small-scale turbulent mass motions. We show that the spatial location of, and timing of the peak in, turbulent kinetic energy together provide persuasive evidence that turbulent energy may play a key role in the transfer of energy in solar flares. Although the kinetic energy of turbulent motions accounts, at any given time, for only ∼ (0.5 -1)% of the energy released, its relatively rapid (∼1 -10 s) energization and dissipation causes the associated throughput of energy (i.e., power) to rival that of major components of the released energy in solar flares, and thus presumably in other astrophysical acceleration sites.

During a solar flare, up to 10 32 ergs of energy stored in magnetic fields in the the solar corona is converted into the energy of accelerated particles, bulk flows, and heating [START_REF] Tandberg-Hanssen | The physics of solar flares[END_REF]. Observations [START_REF] Benz | [END_REF][3][4][5][6] lend considerable support to a scenario in which a significant fraction of the released energy is channelled into accelerated electrons which, guided by the surrounding magnetic field, propagate downward toward the solar surface, producing bremsstrahlung hard X-ray (HXR) emission in collisions with ambient ions along their path [7] and heating the surrounding atmosphere through collisions with ambient electrons [START_REF] Tandberg-Hanssen | The physics of solar flares[END_REF]. This heating of the lower (chromospheric) layers of the solar atmosphere in turn leads to enhanced radiation in Extreme Ultraviolet (EUV) and optical wavelengths and, as a result of the associated increase in gas pressure, to an upward motion of material into the corona [8].

Plasma motions (both inflows and outflows) driven by the primary magnetic reconnection process [START_REF] Priest | Magnetic Reconnection[END_REF] are also observed, both spectroscopically [START_REF] Harra | [END_REF] and through reconfiguration of the magnetic field geometry [11]. The Reynolds number in the solar corona is, as in most astrophysical environments, very large, and accordingly it is expected that these flows will be turbulent [12,13]. Theoretical studies [14] and numerical simulations of magnetic reconnection, on both fluid [15] and kinetic [16,17] scales, have suggested that turbulence can dramatically affect the dynamics of the reconnection process. Further, magnetohydrodynamic (MHD) turbulence has long been conjectured to play a key role in the acceleration of particles during flares [18][START_REF] Ramaty | Particle Acceleration Mechanisms in Astrophysics[END_REF][START_REF] Larosa | [END_REF][21], and numerous models of turbulent (stochastic) acceleration [22][23][24] have been proposed. Together, the above strongly suggests a scenario in which MHD turbulence generated during magnetic reconnection plays a key role in the acceleration of particles; however, to date little firm observational evidence in support of such a scenario has been presented. In this letter, we present multi-faceted observations of an unusually well-observed solar flare that allow an evaluation of the energy content in turbulent plasma motions and hence of the role of such motions in the conversion of magnetic energy to acceleration of fast particles.

A moderately large (GOES class X1.2; see [25]) flare occurred on 2013 May 15 in NOAA solar active region 11748. This flare was observed by several instruments: (1) the Ramaty High Energy Solar Spectroscopic Imager (RHESSI) [26], which produces full-Sun high-spatialresolution soft X-ray (SXR) and hard X-ray (HXR) images with ∼1 keV spectral resolution and several-second time resolution; (2) the Solar Dynamics Observatory (SDO) Helioseismic and Magnetic Imager (HMI) [27], which measures the magnetic field in the lower atmospheric levels; (3) the SDO Atmospheric Imaging Assembly (AIA) [28], which provides high-resolution EUV spatial images; (4) the Hinode EUV Imaging Spectrometer (EIS) [29], which produces EUV spectral line profiles; and (5) the Nobeyama Radioheliograph and Radiopolarimeters [30], which measure the radio-wave radiation produced by mildly relativistic electrons and so provides a diagnosis of the magnetic field strength in the corona.

Fig. 1 shows the general morphology of the flare consistent with the flare reconnection geometry [11]; its nearlimb location allows us to readily ascertain its vertical structure. The flare has a cusp-shaped coronal structure, clearly visible in the AIA 193 Å image which delineates hot flare plasma with temperature ∼ 10 ( 25 keV) footpoints where the accelerated electrons, traveling along the magnetic field lines, impact the relatively dense chromospheric layers of the atmosphere. Using the EUV and SXR images, we estimate the height of the magnetic loop to be ∼1.5 × 10 4 km. Fig. 2 shows the temporal evolution of the flare in radio, SXRs and HXRs. The SXR emission has a spectral shape consistent with bremsstrahlung from a Maxwellian distribution of electrons in a thermal plasma with T ∼ 10 7.5 K, while the HXR emission is produced by bremsstrahlung from higher energy ( 25 keV) nonthermal electrons and typically has a harder, power-law, spectral shape. The peak in HXR emission occurred at ∼01:41 UT, followed at ∼01:45 UT by the peak in SXR emission. The main microwave peaks at 17 and 34 GHz were observed by Nobeyama at ∼01:41 UT, near the time of the peak in HXRs.

Following the methods described in [3,7], which include consideration of the primary bremsstrahlung emission mechanism, instrumental pulse pile-up effects, and the albedo flux resulting from photospheric reflection of primary source HXRs, HXR observations allow us to infer the spatial and spectral distributions of the emitting energetic electrons (Fig. 2). The rate of production, Ṅ , of accelerated electrons above a specified energy E is roughly proportional to the overall intensity of the HXR spectrum above that energy. Because the HXR spectrum, and hence the accelerated electron spectrum that produces it, is quite steep (∝ E -δ with a power-law index δ typically ∼4-6), the total energy in accelerated electrons depends on the shape of the low-energy end of the HXR spectrum. We selected time intervals to avoid instrumental effects such as RHESSI shutter motions and applied the methodology in [31] to produce a range of values of Ṅ consistent with data at each time interval throughout the flare, as shown by the vertical extent of the boxes in Fig. 2.

The emission measure EM=n 2 V and temperature T of the hot thermal SXR-emitting plasma are determined from an isothermal fit [3] to the SXR spectral component; their variations throughout the event are shown in Fig. 2. Using the inferred value of EM and the source volume V estimated from Fig. 1 leads to a estimate of the source density n, which is the lower limit because of the possibility that the emission originates only from a fraction of the observed flare volume, although the estimates [32] indicate that this "filling factor" is consistent with unity.

Broadening of spectral lines in excess of the thermal Doppler width [33] is a signature of turbulent motions, associated either with plasma oscillations [23] or unre-solved bulk plasma motions [8]. To estimate the extent of such turbulence, we use the EIS Fe XXIV spectral line profiles at 192 Å and 255 Å. Since the Fe XXIV 192 Å line represents some 80% of the total intensity in the AIA 193 Å passband, the EIS and AIA images were co-aligned by cross-correlating the EIS 192 Å intensity maps with the AIA 193 Å images. Then, the EIS Fe XXIV 255 Å line profile at each point in the image was fitted using a Gaussian, following the procedure in [START_REF] Harra | [END_REF], which allows for instrumental effects. The characteristic non-thermal broadening velocity v nth is then determined from the extent to which the observed spectral line width exceeds that expected from thermal line broadening [START_REF] Harra | [END_REF]. Fig. 3 shows the spatial distribution of the 255 Å non-thermal line-broadening velocities throughout the source for the time interval starting 01 : 35 : 56 UT (Fig. 2), which corresponds to the interesting epoch just prior to the peak in the HXR light curve. Pixels were excluded where the intensity was either too strong (saturated) or too weak for v nth to be reliably determined. The turbulent velocity tends to be larger by ∼50% near the apex of the magnetic loops and along the outer edge of the arcade [34]. Fig. 2 shows the time variation of v nth , the value averaged over the area A inside the 50% contour of the RHESSI HXR (6 -15 keV) map (Fig. 3).

A similar procedure was used for the Fe XXIV 192 Å line; however, this line was more strongly saturated [34] and hence useful measurements were available only near the start and the end of the flare. Where information from both Fe XXIV 255 Å and 192 Å lines were available (i.e., before 01:36 UT and after 02:03 UT), the inferred values of the average non-thermal broadening velocity v nth agreed within 10% (Fig. 2). Typical values of v nth in this 10 7.2 K plasma were found to be (60 -100) km s -1 . This is somewhat lower than the previously reported (spatially-unresolved) measurements of v nth ≃ 200 km s -1 at higher temperatures [32,33,35,36], suggesting, not surprisingly, that hotter plasma may admit higher turbulent velocities. The total turbulent kinetic energy K ∝ v nth 2 could therefore be larger by a factor of ∼4 than that inferred from the 10 7.2 K lines alone.

The power in non-thermal electrons is given by P = ((δ-1)/(δ-2)) Ṅ E c ; its time history closely matches that of the HXR flux. RHESSI images (Fig. 1) show both the location and the area A of the coronal source, deduced from the 50% intensity contour in the (6-15) keV map, allowing an estimate of the source volume V = A 3/2 ≃ 2 × 10 27 cm 3 , and hence [4] the thermal plasma energy in the coronal source

U th = 3kT √ EM • V , where k is Boltzmann's constant.
The turbulent kinetic energy K at each EIS raster time is calculated using K = (3/2)m i v nth 2 n p V , where m i = 1.3 m p is the mean ion mass for solar coronal abundances [37] and n = EM/V is the number density. Fig. 4 shows the ratio of the turbulent kinetic energy K (ergs) to the instantaneous thermal energy content U th (ergs); the temporal behavior of K is similar to that FIG. 4: Flare energetics. Left panels, top to bottom: power P (erg s -1 ) in non-thermal electrons above the low-energy cutoff Ec, density n (cm -3 ) of the SXR-emitting plasma, thermal energy content U th (erg) of the SXR-emitting plasma. Right panels, top to bottom: bulk kinetic energy K (erg), ratio K/P (s), ratio K/U th (dimensionless).

of U th , with K some two to three orders of magnitude smaller, varying between 10 27 and 10 28 ergs. Since both K and U th ∝ nV ∼ √ EM • V , the effect of a volumetric filling factor less than unity is to reduce them both somewhat, but the ratio K/U th is preserved. While K attains its peak value around the same time (∼01:25 UT) as the SXR flux, it notably has a value equal to some 20% of its peak value as early as 01:40 UT, well before the peak in the HXR flux. Similar behavior is also seen in MHD simulations [38]. Fig. 4 also shows the ratio of the turbulent kinetic energy K (ergs) to the power P (erg s -1 ) in energetic electrons; the ratio K/P (which is a measure of the time it takes a power P to energize/deplete a reservoir of energy K) has a relatively steady value of order 1-10 s.

To estimate the available energy in the magnetic field, we used two independent techniques: (i) microwave spectral data from NoRH and NoRP, and (ii) extrapolated HMI line-of-sight photospheric magnetograms. The microwave spectra were fitted assuming isotropic electrons with a power-law energy spectrum as determined from the RHESSI HXR spectrum. Using fast gyrosynchrotron codes [39], we reproduced the observed NoRP microwave fluxes at three frequencies: 17 GHz and 34 GHz (in the optically thin range) and 9.4 GHz (near the spectral peak); the best-fit spectra corresponded to an av-erage magnetic field strength B ≃ (300 -400) G. The coronal magnetic field strength was also estimated from potential-field (minimum magnetic field strength) extrapolation of the observed HMI line-of-sight photospheric magnetograms, giving B ≃ 300 G at heights ∼ 1.5×10 4 km, where the bulk of the radio emission is observed. These mutually consistent values of the magnetic field strength B correspond to a total magnetic energy (B 2 /8π) V ≃ (7 -12) × 10 30 erg. Following [4], we estimate that the magnetic energy available for dissipation (i.e., the excess over the potential field energy) is 30% of the total magnetic field energy, or ∼ 2 × 10 30 erg.

An enduring challenge in flare physics relates to how such a large fraction of the stored magnetic energy is converted to energy in accelerated particles. In relation to particle acceleration, there are, broadly-speaking, two representations of turbulence (stochasticity): "wave turbulence" [START_REF] Larosa | [END_REF]40,41] and a "stochastic-ensemble-ofcurrent-sheets" [38,[42][43][44]. These two concepts are not necessarily unrelated, since there is a tendency for MHD turbulence to form current sheets [45]. In both scenarios, energy produced at large scales systematically cascades to smaller and smaller scales, where the energy is eventually dissipated to produce heating and acceleration of non-thermal particles. The rate of energy release at large scales and the rate of subsequent energy transfer to smaller, dissipative, scales together determine the rate at which particle acceleration can occur.

In light of this discussion, two aspects of the turbulent energy content K inferred herein are significant. First, the turbulent energy is observed (Figs. 1 and3) to be spatially concentrated in the coronal part of the magnetic loop below the observed cusp-like structure, where the primary energy release is believed to occur. Second, its energy content K grows to a significant level well before the peak in HXR intensity, i.e., before the maximum rate of electron acceleration (Fig. 4). Together, these features lead us to propose that turbulence constitutes a viable channel for the conduit of cascading energy. Although the instantaneous turbulent energy content K is only a percent or so of the available magnetic energy (and of the thermal energy U th in the SXR-emitting plasma), the transfer of energy out of the turbulent energy reservoir could be sufficiently rapid for the associated power to rival that associated with dissipation of the turbulence and the acceleration of non-thermal particles. The ratio of K/P (Fig. 4) shows that for such a scenario to be viable the turbulent energy must be dissipated (and replenished) on a timescale ∼1-10 s. Such a timescale is consistent not only with observed fluctuations in the time profile of the HXR emission in the event studied here, but also with many previous studies [START_REF] Tandberg-Hanssen | The physics of solar flares[END_REF]3].

It is well known [46] that dissipation of anisotropic Alfvén MHD turbulence occurs on a timescale L ⊥ / v nth , where L ⊥ is the characteristic scale associated with variations δB perpendicular to the guiding magnetic field. The "side-on" geometry of this particular flare on the sky (Fig. 1) suggests that the observed line-of-sight velocity fluctuations v nth correspond to motions perpendicular to the guiding field. And although L ⊥ is not directly observable, the dissipation timescale can nevertheless be estimated as follows [46]. The energy density associated with a turbulence-perturbed magnetic field δB is U B ≃ (δB) 2 /8π. Equating this to the turbulent energy content K = (1/2)n m v 2 nth , we obtain v 2 nth ≃ (δB) 2 /4πn m. Since the Alfvén speed V A = B 2 /4πn m, it follows that v nth /V A ≃ δB/B ≃ L ⊥ /L , where L is the longitudinal extent of the turbulence region. Thus the dissipation timescale L ⊥ / v nth is approximately the same as the Alfvén crossing time L /V A , a quantity that is readily ascertainable from observations. Using the inferred values of B and n gives V A ≃ 2 × 10 3 km s -1 for this flare, a typical value for the flaring corona [3]. Thus we expect dissipation of turbulent energy to occur on a timescale L /V A ≃ 5 s, a value consistent both with the inferred value of K/P (Fig. 4) and with the timescales typically associated with the acceleration of electrons by MHD wave turbulence [47]. In the stochastic current sheet models, the ratio U B /K could be different from one, but the observed K could be used to test these models.

In summary, the suite of observations presented herein demonstrate the presence, in the acceleration region, of a significant energy reservoir in turbulent plasma motions which correlates well in time with the acceleration of HXR-producing electrons. An instantaneous energy content ∼10 28 ergs, produced and dissipated on a timescale of a few seconds, transfers a steady-state power ∼ (0.1 -1) × 10 28 erg s -1 , rivalling the power in accelerated non-thermal particles. These observations not only enable quantitative testing of turbulence acceleration models; they lend considerable credence to the idea that turbulence acts as a crucial intermediary in the transfer of energy from reconnecting magnetic fields to accelerated particles during solar flares, and therefore presumably in other astrophysical particle acceleration sites.

FIG. 1 :

 1 FIG. 1: Morphology of the flare. Top: SDO/AIA 193 Å image (background); RHESSI X-ray contours at 50% of peak value for 6 -15 keV (red) and 25 -50 keV (blue) energy ranges, EIS Fe XXIV (255 Å) intensity map (white contours at 30% and 75% of peak value), and Nobeyama 34 GHz radio emission (green contours at 30% and 75% of peak value). Bottom: Cartoon showing the different flare elements and the cooling post-flare magnetic loops.

FIG. 2 :

 2 FIG. 2: Temporal evolution of the 2013 May 15 solar flare parameters. Top to bottom: RHESSI X-ray and Nobeyama radio lightcurves, acceleration rate of non-thermal electrons Ṅ (from RHESSI HXR data), plasma temperature T (from RHESSI SXR data), and the non-thermal broadening velocity v nth (from Hinode/EIS) averaged over the area within the 50% (6 -15) keV contour shown in Fig. 1. The grey dotted vertical lines show the beginning and end of each EIS raster time, and the vertical range of each box indicates the uncertainty in the quantity.

FIG. 3 :

 3 FIG. 3: Non-thermal velocity broadening map for the time interval 01:35:56 UT (see Fig. 2). Background: SDO/AIA 193 Å. Grey scale material: EIS Fe XXIV (255 Å) non-thermal broadening velocity map. Red contour: 50% of maximum intensity in 6 -15 keV HXR, blue contour: 50% of maximum intensity at 25 -50 keV HXR.

Acknowledgments

The authors are thankful to the anonymous two referees for constructive comments and criticism and S. Sridhar for useful discussions. Hinode is a Japanese mission developed and launched by ISAS/JAXA, with NAOJ as a domestic partner and NASA and UKSA as international partners. It is operated by these agencies in cooperation with ESA and NSC (Norway). EPK, NLSJ, NHB were supported by a STFC consolidated grant ST/L000741/1. JEP was supported by CONACyT postdoctoral grant 207887. AGE was supported by grant NNX10AT78G from NASA's Goddard Space Flight Center. AAK was supported in part by the RFBR grant 15-02-03717.