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

Eric Cancès, Paul Cazeaux, Mitchell Luskin

► To cite this version:

Eric Cancès, Paul Cazeaux, Mitchell Luskin. Generalized Kubo Formulas for the Transport Properties of Incommensurate 2D Atomic Heterostructures. Journal of Mathematical Physics, 2017, 58 (6), pp.1-29. 10.1063/1.4984041 . hal-01403588v2

HAL Id: hal-01403588

https://inria.hal.science/hal-01403588v2

Submitted on 18 Nov 2017

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Generalized Kubo Formulas for the Transport Properties of Incommensurate 2D Atomic Heterostructures

Eric Cancès

CERMICS, Ecole des Ponts and Inria Paris

6 & 8 avenue Blaise Pascal

77455 Marne-la-Vallée Cedex 2, France

eric.cances@enpc.fr

Paul Cazeaux

Mitchell Luskin

School of Mathematics

School of Mathematics

University of Minnesota

University of Minnesota

Minneapolis, Minnesota 55455

Minneapolis, Minnesota 55455

pcazeaux@umn.edu

luskin@umn.edu

May 19, 2017

## Abstract

We give an exact formulation for the transport coefficients of incommensurate two-dimensional atomic multilayer systems in the tight-binding approximation. This formulation is based upon the $C^*$ algebra framework introduced by Bellissard and collaborators [2, 3] to study aperiodic solids (disordered crystals, quasicrystals, and amorphous materials), notably in the presence of magnetic fields (quantum Hall effect). We also present numerical approximations and test our methods on a one-dimensional incommensurate bilayer system.

## 1 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 [11]. 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 [23]. Although this approach might provide a good approximation in many cases, the error is generally uncontrolled (see [7] 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 [24]. 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 [13].

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 [2, 3] 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.3.

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 tight-binding model that we consider, we obtain a compact parametrization of the tight-binding Hamiltonian by using environment-dependent site and hopping functions [10]. 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 effects of incommensurability in coupled multilayered systems. We chose to introduce here a simple discretization based on periodic supercells and on the Kernel Polynomial Method [26]. 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 different 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 [6, 16].

Previous research on the development of numerical methods to approximate transport properties within the $C^*$-algebra formulation has been done by Prodan [19] for the effects 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 [16]. We note that the density of states is defined in [16] as a thermodynamic limit, while the density of states and transport properties are given explicit expressions in the
<figure><img src="image_2.png" /><figcaption>Figure 1: Example of a heterostructure: stack of five monolayers</figcaption></figure>

$C^*$-algebra approach, which directly provide the values of the quantities of interest in the thermodynamic limit.

# 2 Geometries of multi-layered systems

## 2.1 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 different 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 $(d+1)$-dimensional space. Note that we choose $d$ as the natural dimension of the structure, since the $(d+1)^{\text{st}}$ dimension plays a very particular and limited role, especially in tight-binding 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 $(d+1)$-dimensional system of $p$ parallel $d$-dimensional periodic atomic layers denoted $\mathcal{L}_j \subset \mathbb{R}^{d+1}$, $j = 1 \dots p$. We denote by

- $(\mathbf{e}_1, \dots, \mathbf{e}_{d+1})$ an orthonormal basis of the physical space such that each layer is perpendicular to $\mathbf{e}_{d+1}$; from now on we identify the physical space with $\mathbb{R}^{d+1}$ using the cartesian coordinates $\mathbf{x} = (x_1, \dots, x_{d+1})^T$ associated with this basis;

- $h_j$ the $(d+1)^{\text{st}}$ coordinate of the center of layer $j$. Without loss of generality, we can assume that $0 = h_1 < h_2 < \dots < h_p$;

- $\mathcal{R}_j$ the $d$-dimensional periodic lattice of layer $j$;
- $\mathbf{E}_j$ the matrix in $\mathbb{R}^{d \times d}$ whose columns form a basis generating the lattice $\mathcal{R}_j$:

$$ \mathcal{R}_j := \mathbf{E}_j \mathbb{Z}^d \subset \mathbb{R}^d; \qquad (1) $$

- $\Gamma_j := \mathbb{R}^d/\mathcal{R}_j$ the quotient of $\mathbb{R}^d$ by the discrete lattice $\mathcal{R}_j$, which has the topology of a $d$-dimensional torus and can be canonically identified with the periodic unit cell $\hat{\Gamma}_j := \mathbf{E}_j[-1/2, 1/2)^d$ of layer $j$;

- $m_j$ the motif of layer $j$, i.e., the measure on $\mathbb{R}^{d+1}$ supported in $\hat{\Gamma}_j \times \mathbb{R}$ and representing the nuclear distribution $\rho_j^{\text{nuc}}$ of layer $j$. More precisely, $m_j$ is a finite sum of positively weighted Dirac measures of the form

$$ m_j = \sum_{k=1}^{M_j} z_k^{(j)} \delta_{\mathbf{x}_k^{(j)}}, $$

where $M_j$ is the number of nuclei per unit cell in layer $j$, $z_1^{(j)}, \dots, z_{M_j}^{(j)}$ the atomic charges of these nuclei, and $\mathbf{x}_k^{(j)} \in \hat{\Gamma}_j \times \mathbb{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 $\rho_j^{\text{nuc}}$ of layer $j$ is then given by

$$ \rho_j^{\text{nuc}} = \sum_{\mathbf{n}_j \in \mathcal{R}_j} m_j \left( \cdot - \left( \gamma_j + \mathbf{n}_j + h_j \mathbf{e}_{d+1} \right) \right) = \sum_{\mathbf{n}_j \in \mathcal{R}_j} \sum_{k=1}^{M_j} z_k^{(j)} \delta_{\gamma_j + \mathbf{n}_j + \mathbf{x}_k^{(j)} + h_j \mathbf{e}_{d+1}}, $$

where $\gamma_j \in \hat{\Gamma}_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 $\mathbf{y}$ of $\mathbb{R}^d$ and its canonical embedding $(\mathbf{y}^T, 0)^T$ in $\mathbb{R}^{d+1}$.

**Remark 2.1.** Note that, by virtue of the identification $\Gamma_j \equiv \hat{\Gamma}_j$, an equivalence class $\gamma_j \in \Gamma_j$ can be seen either as a discrete set of points in $\mathbb{R}^d$, or as one point of the periodic unit cell $\hat{\Gamma}_j$. Since both viewpoints are employed here, we denote these two usages differently to avoid confusion:

- $\gamma_j$ will be used to denote a single point of the periodic unit cell $\hat{\Gamma}_j$,

- $\gamma_j + \mathcal{R}_j$ will be used to denote the corresponding set of points in $\mathbb{R}^d$. In particular, the sum of the values of a function $f$ over the lattice sites of layer $j$ will be denoted as $\sum_{\mathbf{p}_j \in \gamma_j + \mathcal{R}_j} f(\mathbf{p}_j)$.

## 2.2 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 [9]. 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 [2], is simply the periodic unit cell. Indeed, choosing the origin at points
<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of the process for labeling the local configurations of heterostructures. The vertices of the parallelograms are sites of the shifted lattices $\gamma_j + \mathcal{R}_j$.</figcaption></figure>

which differ 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 $\mathbf{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 $\mathbb{R}^d$ on $\mathbb{R}^{d+1}$ by translations $\mathbf{T}$ that are parallel to the layers:

$$ \text{For } \mathbf{a} \in \mathbb{R}^d, \quad \mathbf{T}_{\mathbf{a}} : \begin{cases} \mathbb{R}^{d+1} \to \mathbb{R}^{d+1}, \\ \mathbf{x} \mapsto (x_1 + a_1, \dots, x_d + a_d, x_{d+1}). \end{cases} \qquad (2) $$

**Remark 2.2.** Note that there is no translational symmetry in the perpendicular $(d+1)^{st}$-dimension, and therefore we need not include this direction in the continuous translation group. The situation is different 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 $\mathfrak{M}(\mathbb{R}^{d+1})$ [2], on which the translation group $\mathbb{R}^d$ acts naturally [17]. Then the hull is the dynamical system $(\Omega, \mathbb{R}^d, \mathbf{T})$, where $\Omega$ is the closure of the orbit of the nuclear charge distribution measure on $\mathbb{R}^{d+1}$ generated by the atoms of all $p$ layers under the action of $\mathbb{R}^d$ through $\mathbf{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 efficiently encoded. Since each individual layer $\mathcal{L}_j$, $j = 1, \dots, p$ is periodic, the set of all possible configurations for layer $j$ is in one-to-one correspondence with $\Gamma_j$. Indeed, an element $\gamma_j \in \Gamma_j$ is the projection of the $j$-th layer periodic lattice on the horizontal plane, i.e., an equivalence class modulo $\mathcal{R}_j$, as seen in Figure 2b. The overall configuration can
<figure><img src="image_4.png" /><figcaption>Figure 3: Schematic representation of a configuration $\omega = (\gamma_1, \gamma_2, \gamma_3, \gamma_4) \in \Omega$ for $d = 1$ and $p = 4$. The vector $\gamma_j$ on the figure corresponds to the unique representative of the equivalence class $\gamma_j$ in the periodic unit cell $\hat{\Gamma}_j$ and $\gamma_1 + \mathcal{R}_1$ is the set of the lattice points of the first layer in the configuration $\omega$.</figcaption></figure>

thus be parameterized as an element of

$$ \Omega = \Gamma_1 \times \cdots \times \Gamma_p. \qquad (3) $$

The hull $\Omega$ therefore has the topology of a $(dp)$-dimensional torus. For any given configuration $\omega = (\gamma_1, \dots, \gamma_p)$, the nuclear charge distribution is

$$ \rho_{\omega}^{\text{nuc}} = \sum_{j=1}^{p} \left( \sum_{\mathbf{p}_j \in \gamma_j + \mathcal{R}_j} \mathbf{T}_{\mathbf{p}_j} m_j(\cdot - h_j \mathbf{e}_{d+1}) \right). $$

The support of the measure $\rho_{\omega}^{\text{nuc}}$ is the discrete set:

$$ \begin{cases} \mathcal{L}^{\omega} := \bigcup_{j=1}^{p} \mathcal{L}_{j}^{\omega} \subset \mathbb{R}^{d+1} & \text{(set of all atomic positions),} \\ \text{with } \mathcal{L}_{j}^{\omega} := \text{Supp}(m_j) + \gamma_j + \mathcal{R}_j + h_j \mathbf{e}_{d+1} & \text{(atomic positions in layer j).} \end{cases} \qquad (4) $$

There exists a natural action of the additive $\mathbb{R}^d$ group on $\Omega$, corresponding to an horizontal translation of the whole system in the previous parameterization, defined by

$$ \mathbf{T}_{\mathbf{a}}(\gamma_1, \dots, \gamma_p) = (\gamma_1 + \mathbf{a}, \dots, \gamma_p + \mathbf{a}). \qquad (5) $$

For convenience, we use the same notation $\mathbf{T}$ to denote the action of the translation group $\mathbb{R}^d$ on $\mathbb{R}^{d+1}$, formula (2), and on $\Omega$, formula (5).

Finally, the dynamical system $(\Omega, \mathbb{R}^d, \mathbf{T})$ is equipped with a probability measure $\mathbb{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 [18] lattice $\mathcal{R}$ by $\mathcal{R}^*$.
**Definition 2.3.** The collection of cocompact lattices $\mathcal{R}_1, \dots, \mathcal{R}_p$ of $\mathbb{R}^d$ is called incommensurate if we have for any p-tuple $(\mathbf{k}_1, \dots, \mathbf{k}_p) \in \mathcal{R}_1^* \times \dots \times \mathcal{R}_p^*$,

$$ \sum_{j=1}^{p} \mathbf{k}_j = \mathbf{0} \quad \Leftrightarrow \quad \mathbf{k}_j = \mathbf{0} \quad \forall j = 1, \dots, p. \qquad (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 $\mathcal{R}_1, \dots, \mathcal{R}_p$ of $\mathbb{R}^d$ be cocompact lattices of $\mathbb{R}^d$, $\Gamma_j = \mathbb{R}^d/\mathcal{R}_j$, and $\Omega = \Gamma_1 \times \dots \times \Gamma_p$, endowed with the uniform probability measure $\mathbb{P}$. Then,

1. $\mathbb{P}$ is invariant by the translation group $\mathbb{R}^d$;

2. the dynamical system $(\Omega, \mathbb{R}^d, \mathcal{T}, \mathbb{P})$ is uniquely ergodic if and only if the lattices $\mathcal{R}_1, \dots, \mathcal{R}_p$ are incommensurate. In this case, we have the Birkhoff property: for any $f \in C(\Omega)$ and $\omega \in \Omega$,

$$ \lim_{r \to \infty} \frac{1}{|B_r|} \int_{B_r} f(\mathcal{T}_{-a}\omega) \mathrm{d}\mathbf{a} = \int_{\Omega} f \mathrm{d}\mathbb{P}, \qquad (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 $\mathcal{R}$ be a lattice in $\mathbb{R}^d$ and $\mathbf{k} \in \mathbb{R}^d$. Then:

$$ \lim_{r \to \infty} \frac{1}{\#(\mathcal{R} \cap B_r)} \sum_{\mathbf{n} \in \mathcal{R} \cap B_r} e^{2i\pi \mathbf{k} \cdot \mathbf{n}} = 0 \quad \Leftrightarrow \quad \mathbf{k} \notin \mathcal{R}^*. \qquad (8) $$

It is clear that if $\mathbf{k} \in \mathcal{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 $(\mathbf{c}_1, \dots, \mathbf{c}_d)$ be a basis of $\mathcal{R}$ and $(\mathbf{c}_1^*, \dots, \mathbf{c}_d^*)$ the associated dual basis of $\mathcal{R}^*$ $(\mathbf{c}_j \cdot \mathbf{c}_k^* = \delta_{jk})$. We assume now that $\mathbf{k} \notin \mathcal{R}^*$, and expanding $\mathbf{k} = k_1 \mathbf{c}_1^* + \dots + k_d \mathbf{c}_d^*$, we suppose that $k_1 \notin \mathbb{Z}$ without loss of generality.

Let $\mathcal{E}_r^\perp$ be the projection of $\mathcal{R} \cap B_r$ onto the last $(d-1)$ lattice coordinates:

$$ \mathcal{E}_r^\perp = \left\{ (n_2, \dots, n_d) \in \mathbb{Z}^{d-1} \mid \exists n_1 \in \mathbb{Z}, \sum_{j=1}^d n_j \mathbf{c}_j \in \mathcal{R} \cap B_r \right\}. $$

There exists for any $(d-1)$-tuple $\tilde{\mathbf{n}} = (n_2, \dots, n_d)$ in $\mathcal{E}_r^\perp$ two integers $N_r^\pm(\tilde{\mathbf{n}})$ such that for any $n_1 \in \mathbb{Z}$, $\sum_{j=1}^d n_j \mathbf{c}_j$ belongs to $\mathcal{R} \cap B_r$ if and only if $N_r^-(\tilde{\mathbf{n}}) \le n_1 < N_r^+(\tilde{\mathbf{n}})$. We decompose accordingly the sum in (8):

$$ \sum_{\mathbf{n} \in \mathcal{R} \cap B_r} e^{2i\pi \mathbf{k} \cdot \mathbf{n}} = \sum_{\tilde{\mathbf{n}} \in \mathcal{E}_r^\perp} \sum_{n_1=N_r^-(\tilde{\mathbf{n}})}^{N_r^+(\tilde{\mathbf{n}})-1} e^{2i\pi(k_1 n_1 + \tilde{\mathbf{k}} \cdot \tilde{\mathbf{n}})} \\ = \sum_{\tilde{\mathbf{n}} \in \mathcal{E}_r^\perp} e^{2i\pi \tilde{\mathbf{k}} \cdot \tilde{\mathbf{n}}} \left( \frac{e^{2i\pi k_1 N_r^-(\tilde{\mathbf{n}})} - e^{2i\pi k_1 N_r^+(\tilde{\mathbf{n}})}}{1 - e^{2i\pi k_1}} \right). $$
We can bound the number of elements in $(\mathcal{R} \cap B_r)$ from below by $c_1 r^d$ and of $\mathcal{E}_r^\perp$ by $c_2 r^{d-1}$ where $c_1, c_2 > 0$ are two geometrical constants depending only on $\mathcal{R}$. Thus

$$ \left| \frac{1}{\#(\mathcal{R} \cap B_r)} \sum_{\mathbf{n} \in \mathcal{R} \cap B_r} e^{2i\pi \mathbf{k} \cdot \mathbf{n}} \right| \le \left( \frac{2}{|1 - e^{2i\pi k_1}|} \right) \frac{c_2}{c_1 r}. $$

Since $k_1 \notin \mathbb{Z}$, Lemma 2.5 is proved.

**Birkhoff property.** Let us now suppose that the collection of lattices $\mathcal{R}_1, \dots, \mathcal{R}_p$ is incommensurate. Let $f \in C(\Omega)$ and $\omega = (\gamma_1, \dots, \gamma_p) \in \Omega$. Let $\varepsilon > 0$ and $T_\varepsilon$ be a trigonometric polynomial such that $\|f - T_\varepsilon\|_\infty \le \varepsilon$. $T_\varepsilon$ is a finite linear combination of Fourier factors of the form

$$ G_{\mathbf{k}_1, \dots, \mathbf{k}_p} : (\gamma_1, \dots, \gamma_p) \mapsto e^{2i\pi(\mathbf{k}_1 \cdot \gamma_1 + \dots + \mathbf{k}_p \cdot \gamma_p)} \quad \text{with} \quad (\mathbf{k}_1, \dots, \mathbf{k}_p) \in \mathcal{R}_1^* \times \dots \times \mathcal{R}_p^*. $$

If $\mathbf{k}_j = \mathbf{0}$ for all $j = 1, \dots, p$, then clearly

$$ \lim_{r \to \infty} \frac{1}{|B_r|} \int_{B_r} G_{0, \dots, 0}(T_{-\mathbf{a}}\omega) \mathrm{d}\mathbf{a} = 1 = \int_{\Omega} G_{0, \dots, 0} \mathrm{d}\mathbb{P}. $$

Otherwise, we assume without loss of generality that $\mathbf{k}_2 \ne \mathbf{0}$, and we approximate the ball $B_r$ for $r > 0$ by a union of translated unit cells,

$$ B_r^1 = \bigcup_{\mathbf{n}_1 \in \mathcal{R}_1 \cap B_r} \left( \hat{\Gamma}_1 + \mathbf{n}_1 \right) \quad \text{where} \quad \hat{\Gamma}_1 = \mathbf{E}_1[-1/2, 1/2)^d. $$

Averaging over this approximate ball, we find:

$$ \frac{1}{|B_r^1|} \int_{B_r^1} G_{\mathbf{k}_1, \dots, \mathbf{k}_p} (T_{-\mathbf{a}} \omega) \mathrm{d}\mathbf{a} = \frac{|\hat{\Gamma}_1|^{-1}}{\#(\mathcal{R}_1 \cap B_r)} \sum_{\mathbf{n}_1 \in \mathcal{R}_1 \cap B_r} \int_{\hat{\Gamma}_1} e^{2i\pi \sum_{j=1}^p \mathbf{k}_j \cdot (\gamma_j - \mathbf{a} - \mathbf{n}_1)} \mathrm{d}\mathbf{a} \\ = \left( |\hat{\Gamma}_1|^{-1} \int_{\hat{\Gamma}_1} e^{2i\pi \sum_{j=1}^p \mathbf{k}_j \cdot (\gamma_j - \mathbf{a})} \mathrm{d}\mathbf{a} \right) \frac{1}{\#(\mathcal{R}_1 \cap B_r)} \sum_{\mathbf{n}_1 \in \mathcal{R}_1 \cap B_r} e^{-2i\pi \sum_{j=2}^p \mathbf{k}_j \cdot \mathbf{n}_1}. $$

We assume that the lattices $\mathcal{R}_1, \dots, \mathcal{R}_p$ are incommensurate, and thus $\sum_{j=2}^p \mathbf{k}_j \in \mathcal{R}_1^*$ if and only if $\mathbf{k}_j = \mathbf{0}$ for all $j = 2, \dots, p$. This is not possible since $\mathbf{k}_2 \ne \mathbf{0}$, and we therefore deduce from Lemma 2.5 that, uniformly in $\omega$,

$$ \lim_{r \to \infty} \frac{1}{|B_r^1|} \int_{B_r^1} G_{\mathbf{k}_1, \dots, \mathbf{k}_p} (T_{-\mathbf{a}} \omega) \mathrm{d}\mathbf{a} = 0. $$

Note that for some $C > 0$ independent of $r$, we have $|B_r \Delta B_r^1| \le Cr^{d-1}$ (where $A \Delta B$ denotes the symmetric difference of the sets A and B). Thus we conclude

$$ \forall (\mathbf{k}_1, \dots, \mathbf{k}_p) \in (\mathcal{R}_1^* \times \dots \times \mathcal{R}_p^*) \setminus \{\mathbf{0}\}, \quad \lim_{r \to \infty} \frac{1}{|B_r|} \int_{B_r} G_{\mathbf{k}_1, \dots, \mathbf{k}_p} (T_{-\mathbf{a}} \omega) \mathrm{d}\mathbf{a} = 0, \text{ uniformly in } \omega. $$

As a result, the trigonometric polynomial $T_\varepsilon$ satisfies for $r > 0$ large enough, uniformly in $\omega$,

$$ \left| \frac{1}{|B_r|} \int_{B_r} T_\varepsilon(T_{-\mathbf{a}}\omega) \mathrm{d}\mathbf{a} - \int_{\Omega} T_\varepsilon \mathrm{d}\mathbb{P} \right| \le \varepsilon. $$
Finally, since $\|f - T_\varepsilon\|_\infty \le \varepsilon$ we obtain the Birkhoff ergodic formula: for $r$ large enough, uniformly in $\omega$,

$$ \left| \frac{1}{|B_r|} \int_{B_r} f(\mathbf{T}_{-\mathbf{a}}\omega) \mathrm{d}\mathbf{a} - \int_{\Omega} f \mathrm{d}\mathbb{P} \right| \le 3\varepsilon. \quad (9) $$

**Ergodicity.** Let $B$ be a measurable subset of $\Omega$ invariant under the action of $\mathbb{R}^d$. Let $\varepsilon > 0$ and $f_\varepsilon$ be a continuous function such that $\|\chi_B - f_\varepsilon\|_{L^1(\Omega)} \le \varepsilon$. By invariance of $B$ under translations, we have for all $\mathbf{a} \in \mathbb{R}^d$,

$$ \|f_\varepsilon \circ \mathbf{T}_\mathbf{a} - f_\varepsilon\|_{L^1(\Omega)} \le 2\varepsilon. $$

We can bound the difference between $f_\varepsilon$ and its Birkhoff means:

$$ \begin{aligned} \left\| \frac{1}{|B_r|} \int_{B_r} f_\varepsilon \circ \mathbf{T}_{-\mathbf{a}} \mathrm{d}\mathbf{a} - f_\varepsilon \right\|_{L^1(\Omega)} &= \left| \int_{\Omega} \left| \frac{1}{|B_r|} \int_{B_r} f_\varepsilon(\mathbf{T}_{-\mathbf{a}}\omega) \mathrm{d}\mathbf{a} - f_\varepsilon(\omega) \right| \mathrm{d}\mathbb{P}(\omega) \right| \\ &\le \frac{1}{|B_r|} \int_{B_r} \int_{\Omega} \left| f_\varepsilon(\mathbf{T}_{-\mathbf{a}}\omega) - f_\varepsilon(\omega) \right| \mathrm{d}\mathbb{P}(\omega) \mathrm{d}\mathbf{a} \le 2\varepsilon, \end{aligned} $$

where we have used the triangle inequality for integrals and Fubini's theorem for non-negative functions. Since the Birkhoff means of $f_\varepsilon$ converge uniformly in $\omega$ for $r \to \infty$ to $\int_{\Omega} f_\varepsilon \mathrm{d}\mathbb{P}$, we deduce that

$$ \left\| \int_{\Omega} f_{\varepsilon} \mathrm{d}\mathbb{P} - f_{\varepsilon} \right\|_{L^1(\Omega)} \le 2\varepsilon. $$

As a consequence, we have by the triangle inequality,

$$ \|\chi_B - \mathbb{P}(B)\|_{L^1(\Omega)} \le \|\chi_B - f_\varepsilon\|_{L^1(\Omega)} + \left\| f_\varepsilon - \int_{\Omega} f_\varepsilon \mathrm{d}\mathbb{P} \right\|_{L^1(\Omega)} + \left\| \int_{\Omega} f_\varepsilon \mathrm{d}\mathbb{P} - \mathbb{P}(B) \right\|_{L^1(\Omega)} \le 4\varepsilon. $$

Since this holds for any $\varepsilon > 0$, we conclude that $\chi_B$ is constant a.s., and therefore $\mathbb{P}(B) \in \{0, 1\}$. Thus $\mathbb{P}$ is ergodic.

**Non-averaging case.** If the lattices $\mathcal{R}_1, \dots, \mathcal{R}_p$ are not incommensurate, there exists a particular nonzero combination $(\mathbf{k}_1, \dots, \mathbf{k}_p) \in \mathcal{R}_1^* \times \dots \times \mathcal{R}_p^*$ such that $\sum_{j=1}^p \mathbf{k}_j = \mathbf{0}$. Then, we have for the corresponding Fourier factor,

$$ G_{\mathbf{k}_1, \dots, \mathbf{k}_p}(\mathbf{T}_{-\mathbf{a}}\omega) = G_{\mathbf{k}_1, \dots, \mathbf{k}_p}(\omega), \quad \forall (\mathbf{a}, \omega) \in \mathbb{R}^d \times \Omega. $$

The function $G_{\mathbf{k}_1, \dots, \mathbf{k}_p}$ is then invariant by the $\mathbb{R}^d$-action and not constant. Therefore the dynamical system $(\Omega, \mathbb{R}^d, \mathbf{T}, \mathbb{P})$ is not ergodic in this case. $\square$

**Remark 2.6.** When the lattices $\mathcal{R}_1, \dots, \mathcal{R}_p$ are commensurate, i.e., $\bigcap_{j=1}^p \mathcal{R}_j$ forms a cocompact superlattice of $\mathbb{R}^d$, then the orbits through the action $\mathbf{T}$ of $\mathbb{R}^d$ on $\Omega$ are no longer dense in $\Omega$, rather they form lower dimensional submanifolds of $\Omega$ 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 $\mathcal{R}_1, \dots, \mathcal{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 $\mathbb{R}^d$ may not be dense in $\Omega$, nor has a simple geometry.
## 2.3 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 $\mathcal{R}_0$-periodic potential, and whose top layer may have defects. We denote by $\mathcal{R}_0$ the periodic lattice of the substrate, and by $\mathcal{R}_1$ and $\mathcal{R}_2$ the periodic lattices of the bottom and top layers respectively, and we assume that $\mathcal{R}_0$, $\mathcal{R}_1$ and $\mathcal{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 $\mathcal{R}_2$ is still appropriate to describe the configurations of the system; the difference 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: $(m_2^{(0)}, m_2^{(1)}, \dots, m_2^{(D)})$, where $m_2^{(0)}$ is the periodic motif of the top layer in the absence of defects, and $m_2^{(1)}, \dots, m_2^{(D)}$ correspond to the different 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_2^{(0)}$ corresponds to the case when no hydrogen atom is adsorbed, $m_2^{(1)}$ (respectively, $m_2^{(2)}$) to the case when only the first (respectively, second) carbon atom has adsorbed a hydrogen atom, and $m_2^{(3)}$ 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 $(p_0, p_1, \dots, p_D)$ (with $p_j > 0$ and $\sum_{j=0}^D p_j = 1$) the probability distribution of the motifs $(m_2^{(0)}, m_2^{(1)}, \dots, m_2^{(D)})$.

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 $\mathbf{r} \in \mathbb{R}^d$, we introduce the decomposition of $\mathbf{r}$ associated with the lattice $\mathcal{R}_2$ defined as

$$ \mathbf{r} = [\mathbf{r}]_2 + \{\mathbf{r}\}_2 \quad \text{with } [\mathbf{r}]_2 \in \mathcal{R}_2 \text{ and } \{\mathbf{r}\}_2 \in \hat{\Gamma}_2 = \mathbf{E}_2[-1/2, 1/2)^d. \quad (10) $$

For this example, the hull is the dynamical system $(\Omega, \mathbb{R}^d, \tau, \mathbb{P})$ where

- $\Omega = \Gamma_0 \times \Gamma_1 \times \Gamma_2 \times \{0, 1, \dots, D\}^{\mathcal{R}_2}$, the configuration $\omega = (\gamma_0, \gamma_1, \gamma_2, \{t_{\mathbf{n}_2}\}_{\mathbf{n}_2 \in \mathcal{R}_2}) \in \Omega$ corresponding to the nuclear distribution

$$ \rho_{\omega}^{\text{nuc}} = \sum_{\mathbf{p}_1 \in \gamma_1 + \mathcal{R}_1} T_{\mathbf{p}_1} m_1 + \sum_{\mathbf{p}_2 \in \gamma_2 + \mathcal{R}_2} T_{\mathbf{p}_2} m_2^{(t_{[\mathbf{p}_2]_2})} (-h \mathbf{e}_{d+1}). $$

Here $t_{[\mathbf{p}_2]_2} \in \{0, 1, \dots, D\}$ is the type of the motif carried by the cell of the top layer centered at $\mathbf{p}_2 + h\mathbf{e}_{d+1}$. Note that $\gamma_0$ does not appear in the definition of $\rho_{\omega}^{\text{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 $\mathcal{R}_0$-periodic potential which modifies the electronic structure of the bilayer system, and this potential depends on $\gamma_0$, the relative position of the substrate with respect to the bilayer system;
- $\tau$ is the action of the group $\mathbb{R}^d$ on $\Omega$ defined by

$$ \forall \mathbf{a} \in \mathbb{R}^d, \quad \forall \omega = (\gamma_0, \gamma_1, \gamma_2, \{t_{\mathbf{n}_2}\}_{\mathbf{n}_2 \in \mathcal{R}_2}), \\ \tau_{\mathbf{a}}(\omega) = (\gamma_0 + \mathbf{a}, \gamma_1 + \mathbf{a}, \gamma_2 + \mathbf{a}, \{t_{\mathbf{n}_2 - \{\gamma_2\}_{2+\mathbf{a}}\}_{\mathbf{n}_2 \in \mathcal{R}_2}\}); $$

- denoting by $\mu$ the probability on the set $\{0, 1, \dots, D\}$ with law $(p_0, p_1, \dots, p_D)$, the ergodic probability $\mathbb{P}$ on $\Omega$ is defined as

$$ \mathbb{P} = \mathbb{P}_0 \otimes \mathbb{P}_1 \otimes \mathbb{P}_2 \otimes \mathbb{P}_{\text{dis}}, $$

where $\mathbb{P}_j$ is the uniform probability distribution on $\Gamma_j$, and $\mathbb{P}_{\text{dis}} = \mu^{\otimes \mathcal{R}_2}$ is the probability distribution on the disorder.

## 3 $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 [2, 3] 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:

- $\Omega_D = \{1, \dots, p\} \times \Omega$;

- $\mathbf{t}$ is the action of the group

$$ \mathbb{G}_D = \mathbb{Z}/p\mathbb{Z} \times \mathbb{R}^d $$

defined by: for all $a = (\alpha, \mathbf{a}) \in \mathbb{G}_D$, and all $(j, \gamma_1, \dots, \gamma_p) \in \Omega_D$,

$$ \mathbf{t}_a(j, \gamma_1, \dots, \gamma_p) = (j - \alpha, \gamma_1 + \mathbf{a}, \dots, \gamma_p + \mathbf{a}), $$

where the $-$ sign in $j - \alpha$ refers to the natural action of $\mathbb{Z}/p\mathbb{Z}$ on $\{1, \dots, p\}$ (if $\alpha = k + p\mathbb{Z}$, $j - \alpha = ((j - k - 1) \mod p) + 1$);

- $\mathbb{P}_D$ is as usual the uniform probability distribution on $\Omega_D$, and it is ergodic if and only if the lattices $\mathcal{R}_1, \dots, \mathcal{R}_p$ are incommensurate.

More explicitly, for non-disordered multilayers, $\Omega_D$ consists of $p$ copies of $\Gamma_1 \times \dots \times \Gamma_p$ indexed by $j \in \{1, \dots, p\}$ 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 [14, 15, 25]. We therefore introduce, for each layer $k$, a finite set $\Xi_k$ of tight-binding orbitals per unit cell. For a given configuration $\omega = (j, \gamma_1, \dots, \gamma_p)$, the tight-binding orbital of type $n \in \Xi_k$ associated with the lattice site of layer $k$ located at point $\mathbf{x} + h_k \mathbf{e}_{d+1} \in \gamma_k + \mathcal{R}_k + h_k \mathbf{e}_{d+1}$ is indexed by the triplet $((k - j) \mod p, \mathbf{x}, n) \in \mathbb{Z}/p\mathbb{Z} \times (\gamma_k + \mathcal{R}_k) \times \Xi_k$. The integer $((k - j) \mod p)$, considered as a element of $\mathbb{Z}/p\mathbb{Z}$, 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 $\omega$, 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 $\Xi^\omega$ of all tight-binding orbitals of the $p$ layers is then indexed by:

$$ \Xi^\omega = \bigcup_{k=1}^{p} \{(k-j) \mod p\} \times (\gamma_k + \mathcal{R}_k) \times \Xi_k. \quad (11) $$

In the special case when there is only one orbital per unit cell, each $\Xi_k$ only contains one element, and can therefore be omitted in the above definition of $\Xi^\omega$.

In the configuration $\omega$, the space of quantum states for the tight-binding model is

$$ \mathfrak{H}_{\omega} = \ell^2(\Xi^{\omega}), \quad (12) $$

and observables such as the Hamiltonian are described as linear operators on $\mathfrak{H}_\omega$.

## 3.1 Abstract setting

Let us first briefly recall in this section the general formalism of groupoid $C^*$-algebras. We refer e.g. to [3] 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

$$ \Xi^{\omega} \equiv \bigcup_{k=1}^{p} \{(k-j) \mod p\} \times (\gamma_k + \mathcal{R}_k) \subset \mathbb{G}_D. \quad (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 $\Omega_D$:

$$ X = \{\omega \in \Omega_D \mid \omega = (j, \gamma_1, \dots, \gamma_p); \gamma_j = \mathbf{0}\}. \quad (14) $$

The idea is that while $\Omega_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 $\Gamma(X)$ defined as follows:

$$ \Gamma(X) = \{(\omega, a) \in X \times \mathbb{G}_D \mid \mathbf{t}_{-a}\omega \in X\}. \quad (15) $$

For $\omega = (j, \gamma_1, \dots, \gamma_p) \in \Omega_D$ and $a = (\alpha, \mathbf{a}) \in G_D$, $(\omega, a) \in \Gamma(X)$ if and only if $\gamma_j = 0$ and $\gamma_{j+\alpha} = \mathbf{a}$. Geometrically, $\Gamma(X)$ 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 $(\omega, a)$ is $\mathbf{A} = \mathbf{a} + (h_{j+\alpha} - h_j)\mathbf{e}_{d+1}$. The groupoid is equipped with a set of three operations, respectively the range $r: \Gamma(X) \to X$, source $s: \Gamma(X) \to X$, and composition $\circ: \Gamma(X) \times \Gamma(X) \to \Gamma(X)$, satisfying:

$$ r(\omega, a) = \omega, \quad s(\omega, a) = \mathbf{t}_{-a}\omega, \quad (\omega, a) \circ (\mathbf{t}_{-a}\omega, b) = (\omega, a+b). \quad (16) $$

The fiber $\Gamma^{(\omega)}$ is defined as $r^{-1}(\{\omega\})$ for any $\omega \in X$.

In a second step, we define the $*$-algebra $\mathcal{A}_0$ of continuous functions with compact support defined on $\Gamma(X)$ with values in $\mathbb{C}$, endowed with the following composition laws
and *-operator:

$$ (\lambda f)(\omega, a) = \lambda f(\omega, a), \\ (f + g)(\omega, a) = f(\omega, a) + g(\omega, a), \\ (f * g)(\omega, a) = \sum_{(\omega,b) \in \Gamma^{(\omega)}} f(\omega, b)g(\mathbf{t}_{-b}\omega, a-b), \qquad (17) \\ f^*(\omega, a) = \overline{f(\mathbf{t}_{-a}\omega, -a)}. $$

The *-algebra $\mathcal{A}_0$ has an identity, denoted by $\mathbf{1}$:

$$ \mathbf{1}(\omega, a) = \delta_{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 (17) do not take their values in $\mathbb{C}$ but $f(\omega, a)$ is in $\mathbb{C}^{N_j \times N_{j+\alpha}}$ when $\omega = (j, \gamma_1, \dots, \gamma_p)$ and $a = (\alpha, \mathbf{a})$, where $N_k = \#(\Xi_k)$ is the total number of tight-binding orbitals carried by the atoms in the unit cell of layer $k$. Products of the form $f(\omega, b)g(\mathbf{t}_{-b}\omega, a-b)$ should then be understood as matrix products.

The *-algebra $\mathcal{A}_0$ can be mapped onto a *-algebra of bounded linear operators acting on the space of quantum states $\mathfrak{H}_\omega = \ell^2(\Xi^\omega)$, see (12), via the representation formula:

$$ \pi_{\omega}(f)\phi(x) = \sum_{(\omega,y) \in \Gamma^{(\omega)}} f(\mathbf{t}_{-x}\omega, y-x) \phi(y), \quad \forall \phi \in \mathfrak{H}_{\omega}, \quad \forall x \in \Xi^{\omega}. \qquad (18) $$

In the above formula, we have implicitly used the fact that the set $\Xi^\omega$ defined by (13) coincides with $\{y \in \mathbb{G}_D | (\omega, y) \in \Gamma^{(\omega)}\}$.

<figure><img src="image_5.png" /><figcaption>Figure 4: Graphical explanation of the representation formula (18) for $d = 1, p = 4$, and for a configuration $\omega = (3, \gamma_1, \gamma_2, \mathbf{0}, \gamma_4) \in \Omega_D$. The points $x = (2-3, \mathbf{x}) = (-1, \mathbf{x})$ and $y = (4-3, \mathbf{y}) = (1, \mathbf{y})$ of $\Xi^\omega$ are combined to form an element $y-x = (2, \mathbf{y}-\mathbf{x})$ such that $(\mathbf{t}_{-x}\omega, y-x) \in \Gamma^{(\mathbf{t}_{-x}\omega)}$.</figcaption></figure>

The following covariance condition holds: for a given $\gamma = (\omega, a) \in \Gamma(X)$,

$$ \pi_{\mathbf{t}_{-a}\omega}(f) = \mathfrak{T}(\gamma)^{-1}\pi_{\omega}(f)\mathfrak{T}(\gamma), \qquad (19) $$

where the translations $\mathfrak{T}(\gamma) : \mathfrak{H}_{\mathbf{t}_{-a}\omega} \to \mathfrak{H}_{\omega}$ are unitary operators defined by

$$ \mathfrak{T}(\gamma)\phi(x) = \phi(\mathbf{t}_{-a}x), \quad \forall \phi \in \mathfrak{H}_{\mathbf{t}_{-a}\omega}, \quad \forall x \in \Xi^{\omega}. \qquad (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 $\mathcal{B}$ be an
antisymmetric 3 × 3 matrix representing the magnetic field. For $x = (l, \mathbf{x})$ and $y = (m, \mathbf{y})$ in $\Xi^\omega$ with $\omega = (j, \gamma_1, \dots, \gamma_p)$, we define the corresponding Peierls phase term as

$$ \Phi_{\mathcal{B}}^{\omega}(x, y) = \pi \sum_{\nu, \mu=1}^{3} \mathcal{B}_{\nu, \mu} X_{\nu} Y_{\mu}, \qquad (21) $$

where $\mathbf{X} = \mathbf{x} + (h_{j+l} - h_j)\mathbf{e}_3$ and $\mathbf{Y} = \mathbf{y} + (h_{j+m} - h_j)\mathbf{e}_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

$$ \begin{aligned} (f *_{\mathcal{B}} g)(\omega, a) &= \sum_{(\omega, b) \in \Gamma^{(\omega)}} f(\omega, b)g(\mathbf{t}_{-b}\omega, a-b) e^{i\Phi_{\mathcal{B}}^{\omega}(a,b)}, \\ \pi_{\omega}(f)\phi(x) &= \sum_{(\omega, y) \in \Gamma^{(\omega)}} f(\mathbf{t}_{-x}\omega, y-x) e^{-i\Phi_{\mathcal{B}}^{\omega}(x,y)}\phi(y), \quad \forall \phi \in \mathfrak{H}_{\omega}, x \in \Xi^{\omega}, \end{aligned} \qquad (22) $$

where we again used the fact that $\Xi^\omega = \{y \in \mathbb{G}_D \mid (\omega, y) \in \Gamma^{(\omega)}\}$. In this case, the following covariance condition holds: for a given $\gamma = (\omega, a) \in \Gamma(X)$,

$$ \pi_{\mathbf{t}_{-a}\omega}(f) = U(\gamma)^{-1}\pi_{\omega}(f)U(\gamma), \qquad (23) $$

where the magnetic translations $U(\gamma) : \mathfrak{H}_{\mathbf{t}_{-a}\omega} \to \mathfrak{H}_{\omega}$ are unitary operators defined by

$$ U(\gamma)\phi(x) = \exp \left( i \int_{[\mathbf{X}-\mathbf{A}, \mathbf{X}]} \mathcal{A}(y) \cdot dy \right) \phi(\mathbf{t}_{-a}x), \quad \forall \phi \in \mathfrak{H}_{\mathbf{t}_{-a}\omega}, \quad \forall x \in \Xi^{\omega}. \qquad (24) $$

Here, $\mathcal{A}$ is a vector potential giving rise to the magnetic field $\mathcal{B}$, and $[\mathbf{X}-\mathbf{A}, \mathbf{X}]$ is the line segment joining $\mathbf{X}-\mathbf{A}$ to $\mathbf{X}$ in $\mathbb{R}^3$.

Note that $\pi_\omega(f)$ is hermitian if $f = f^*$. This representation induces a $C^*$ norm on $\mathcal{A}_0$, defined by:

$$ \|f\| = \sup_{\omega \in X} \|\pi_{\omega}(f)\|, $$

where the norm on the right hand side is the operator norm on $\mathcal{L}(\mathfrak{H}_\omega)$, the space of bounded linear operators on $\mathfrak{H}_\omega$. We then construct in a third and final step the $C^*$-algebra $\mathcal{A}$ as the completion of $\mathcal{A}_0$ for this norm.

**Integro-differential calculus** Integration and derivations can be constructed on the $C^*$-algebra as follows. Let $\mathbb{P}$ be a probability measure on $X$ invariant by $\Gamma(X)$-action. A formal integration is then obtained on $\mathcal{A}_0$ as follows: if $f \in \mathcal{A}_0$,

$$ \mathcal{T}_{\mathbb{P}}(f) = \int_X d\mathbb{P}(\omega) f(\omega, \mathbf{0}). \qquad (25) $$

This trace is positive, $\mathcal{T}_{\mathbb{P}}(f^* * f) \ge 0$. It is faithful ($\mathcal{T}_{\mathbb{P}}(f^* * f) = 0$ iff $f = 0$) when the support of the measure $\mathbb{P}$ is $X$.

**Remark 3.3.** This trace coincides with the trace per unit volume in $\mathbb{R}^d$ of the corresponding operator $\pi_\omega(f)$ thanks to the Birkhoff property (7).
One then defines $L^p$-norms by setting for $1 \le p < \infty$,

$$ \|f\|_p = \mathcal{T}_p \left( (f * f^*)^{p/2} \right)^{1/p} \quad \forall f \in \mathcal{A}_0. \qquad (26) $$

Closure of $\mathcal{A}_0$ with respect to the $L^p$-norms defines Banach spaces $L^p(\mathcal{A}, \mathcal{T}_p)$. Derivations are introduced next by setting

$$ \partial_j f(\omega, x) = i x_j f(\omega, x), \quad f \in \mathcal{A}_0. \qquad (27) $$

Given a multi-index $\alpha = (\alpha_1, \dots, \alpha_d)$ we will use the notation $\partial_\alpha f$ for $\partial_1^{\alpha_1} \dots \partial_d^{\alpha_d}$. The derivations defined above satisfy the fundamental properties of derivation operators: they commute, they are $*$-derivations in the sense that $\partial_j(f^*) = (\partial_j f)^*$, and satisfy the Leibniz rule $\partial_j(f * g) = \partial_j f * g + f * \partial_j g$. Finally one has the operator representation:

$$ \pi_\omega(\partial_j f) = -i [x_j^\omega, \pi_\omega f], \qquad (28) $$

where $x^\omega = (x_1^\omega, \dots, x_d^\omega)$ is the position operator on $\mathcal{L}^\omega$.

**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 [1] for the details. For a given element $f \in \mathcal{A}$, its resolvent set $\rho(f)$ and spectrum $\sigma(f)$ are

$$ \rho(f) = \{z \in \mathbb{C} \text{ s.t. } (z\mathbf{1} - f) \text{ is invertible}\}, \quad \sigma(f) = \mathbb{C}\setminus\rho(f). $$

Furthermore, the resolvent set is open in $\mathbb{C}$, and the resolvent function $z \mapsto (z\mathbf{1} - f)^{-1}$ is an algebra-valued analytic function of $z \in \rho(f)$. An analytical calculus can then be defined on $\mathcal{A}$ as follows. Let $\mathcal{F}_{\sigma(f)}$ by the algebra of $\mathbb{C}$-valued functions which are analytic in an open neighborhood of $\sigma(f)$. Then we define a homomorphism of algebras

$$ \mathcal{F}_{\sigma(f)} \ni \phi \mapsto \frac{1}{2i\pi} \oint_C \phi(z)(z\mathbf{1} - f)^{-1} dz \in \mathcal{A}, \qquad (29) $$

where $C$ is a contour surrounding $\sigma(f)$ in the analyticity domain of $\phi$. The right-hand side of this mapping is independent of the particular choice of contour, and usually noted as $\phi(f)$.

## 3.2 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 $\Omega_D$ of a perfect incommensurate bilayer. We have

$$ \Omega_D = \{1, 2\} \times \Gamma_1 \times \Gamma_2, $$
and

$$X = \{(1, \mathbf{0}, \gamma_2) ; \gamma_2 \in \Gamma_2\} \cup \{(2, \gamma_1, \mathbf{0}) ; \gamma_1 \in \Gamma_1\}.$$

This leads to the identification:

$$X = X_1 \cup X_2 \quad \text{where} \quad \left\{ \begin{array}{l} X_1 = \Gamma_2, \\ X_2 = \Gamma_1. \end{array} \right. \qquad (30)$$

The set $X_1$ (resp. $X_2$) describes the possible configurations from the point of view of a given lattice site of layer $\mathcal{L}_1$ (resp. $\mathcal{L}_2$). In a configuration $\gamma_2 \in X_1$, which corresponds to $\omega = (1, \mathbf{0}, \gamma_2)$, the lattice sites of $\mathcal{L}_1$ are located at $\mathcal{R}_1$ while the lattice sites of $\mathcal{L}_2$ are located at $\gamma_2 + \mathcal{R}_2 + h\mathbf{e}_{d+1}$, where h is the distance between the two layers. In a configuration $\gamma_1 \in X_2$, which corresponds to $\omega = (2, \gamma_1, \mathbf{0})$, the lattice sites of $\mathcal{L}_2$ are located at $\mathcal{R}_2$ while the lattice sites of $\mathcal{L}_1$ are located at $\gamma_1 + \mathcal{R}_1 - h\mathbf{e}_{d+1}$.

The decomposition (30) of the transversal allows us to identify a block decomposition of the groupoid:

$$\Gamma(X) = \vec{\Gamma}_{11} \cup \vec{\Gamma}_{12} \cup \vec{\Gamma}_{21} \cup \vec{\Gamma}_{22} = \vec{\Gamma}(X_1, X_2), \qquad (31)$$

where the set of arrows $\vec{\Gamma}_{11}$ and $\vec{\Gamma}_{22}$ include all intralayer jumps:

$$\left\{ \begin{array}{l} \vec{\Gamma}_{11} = \{(\gamma_2, \mathbf{m}) ; \gamma_2 \in X_1, \mathbf{m} \in \mathcal{R}_1\}, \\ \vec{\Gamma}_{22} = \{(\gamma_1, \mathbf{n}) ; \gamma_1 \in X_2, \mathbf{n} \in \mathcal{R}_2\}, \end{array} \right.$$

while $\vec{\Gamma}_{21}$ describes jumps from the first to the second layer, and $\vec{\Gamma}_{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 $(\gamma_2, \mathbf{q})$ in $\vec{\Gamma}_{12}$. Indeed, $\gamma_2 \in X_1$ can be deduced from $\mathbf{q} \in \gamma_2 + \mathcal{R}_2$ as its equivalence class modulo $\mathcal{R}_2$. We will thus denote these arrows solely as arrow vectors $\vec{\mathbf{q}}$, observing that these can take their values in all of $\mathbb{R}^d$:

$$\left\{ \begin{array}{l} \vec{\Gamma}_{12} = \{\vec{\mathbf{q}}; \mathbf{q} \in \mathbb{R}^d\}, \\ \vec{\Gamma}_{21} = \{\vec{\mathbf{p}}; \mathbf{p} \in \mathbb{R}^d\}. \end{array} \right. \qquad (32)$$

This allows us to make explicit the three groupoid operations, namely the range map

$$r : \vec{\Gamma} \to \vec{\Gamma}, \quad \left\{ \begin{array}{l} \vec{\Gamma}_{11} \to X_1, \quad (\gamma_2, \mathbf{m}) \mapsto \gamma_2, \\ \vec{\Gamma}_{12} \to X_1, \quad \vec{\mathbf{q}} \mapsto \vec{\mathbf{q}}^{\mathcal{R}_2}, \\ \vec{\Gamma}_{21} \to X_2, \quad \vec{\mathbf{p}} \mapsto \vec{\mathbf{p}}^{\mathcal{R}_1}, \\ \vec{\Gamma}_{22} \to X_2, \quad (\gamma_1, \mathbf{n}) \mapsto \gamma_1, \end{array} \right. \qquad (33)$$

the source map:

$$s : \vec{\Gamma} \to \vec{\Gamma}, \quad \left\{ \begin{array}{l} \vec{\Gamma}_{11} \to X_1, \quad (\gamma_2, \mathbf{m}) \mapsto \gamma_2 - \overline{\mathbf{m}}^{\mathcal{R}_2}, \\ \vec{\Gamma}_{12} \to X_2, \quad \vec{\mathbf{q}} \mapsto -\vec{\mathbf{q}}^{\mathcal{R}_1}, \\ \vec{\Gamma}_{21} \to X_1, \quad \vec{\mathbf{p}} \mapsto -\vec{\mathbf{p}}^{\mathcal{R}_2}, \\ \vec{\Gamma}_{22} \to X_2, \quad (\gamma_1, \mathbf{n}) \mapsto \gamma_1 - \overline{\mathbf{n}}^{\mathcal{R}_1}, \end{array} \right. \qquad (34)$$
and the composition map:

$$
\left\{
\begin{array}{llll}
\text{For } (\gamma_2, \mathbf{m}) \in \bar{\Gamma}_{11}, & \mathbf{m}' \in \mathcal{R}_1 : (\gamma_2, \mathbf{m}) \circ (\gamma_2 - \bar{\mathbf{m}}^{\mathcal{R}_2}, \mathbf{m}') & = (\gamma_2, \mathbf{m} + \mathbf{m}') & \in \bar{\Gamma}_{11}, \\
\text{for } \vec{\mathbf{q}} \in \bar{\Gamma}_{12}, & \mathbf{p} + \mathbf{q} \in \mathcal{R}_1 : \vec{\mathbf{q}} \circ \vec{\mathbf{p}} & = (\vec{\mathbf{q}}^{\mathcal{R}_2}, \mathbf{q} + \mathbf{p}) & \in \bar{\Gamma}_{11}, \\
\text{for } \vec{\mathbf{p}} \in \bar{\Gamma}_{21}, & \mathbf{q} + \mathbf{p} \in \mathcal{R}_2 : \vec{\mathbf{p}} \circ \vec{\mathbf{q}} & = (\vec{\mathbf{p}}^{\mathcal{R}_1}, \mathbf{p} + \mathbf{q}) & \in \bar{\Gamma}_{22}, \\
\text{for } (\gamma_1, \mathbf{n}) \in \bar{\Gamma}_{22}, & \mathbf{n}' \in \mathcal{R}_2 : (\gamma_1, \mathbf{n}) \circ (\gamma_1 - \bar{\mathbf{n}}^{\mathcal{R}_1}, \mathbf{n}') & = (\gamma_1, \mathbf{n} + \mathbf{n}') & \in \bar{\Gamma}_{22}, \\
\text{for } \vec{\mathbf{p}} \in \bar{\Gamma}_{21}, & \mathbf{m} \in \mathcal{R}_1 : \vec{\mathbf{p}} \circ (-\vec{\mathbf{p}}^{\mathcal{R}_2}, \mathbf{m}) & = \overrightarrow{\mathbf{p} + \mathbf{m}} & \in \bar{\Gamma}_{21}, \\
\text{for } (\gamma_1, \mathbf{n}) \in \bar{\Gamma}_{22}, & \mathbf{p} + \mathbf{n} \in \gamma_1 + \mathcal{R}_1 : (\gamma_1, \mathbf{n}) \circ \vec{\mathbf{p}} & = \overrightarrow{\mathbf{n} + \vec{\mathbf{p}}} & \in \bar{\Gamma}_{21}, \\
\text{for } (\gamma_2, \mathbf{n}) \in \bar{\Gamma}_{11}, & \mathbf{q} + \mathbf{n} \in \gamma_2 + \mathcal{R}_2 : (\gamma_2, \mathbf{m}) \circ \vec{\mathbf{q}} & = \overrightarrow{\mathbf{m} + \vec{\mathbf{q}}} & \in \bar{\Gamma}_{12}, \\
\text{for } \vec{\mathbf{q}} \in \bar{\Gamma}_{12}, & \mathbf{n} \in \mathcal{R}_2 : \vec{\mathbf{q}} \circ (-\vec{\mathbf{q}}^{\mathcal{R}_1}, \mathbf{n}) & = \overrightarrow{\mathbf{q} + \vec{\mathbf{n}}} & \in \bar{\Gamma}_{12},
\end{array}
\right.
\tag{35}
$$

where we have introduced the notation $\bar{\mathbf{r}}^{\mathcal{R}_j} = \mathbf{r} + \mathcal{R}_j$, $j = 1 \dots p$ to denote the equivalence classes in $\Gamma_j$ of a given point $\mathbf{r} \in \mathbb{R}^d$. Finally, the fiber $\Gamma^{(\omega)}$ can be determined:

$$
\vec{\Gamma}(\gamma_2) = \Gamma^{(1,0,\gamma_2)} \equiv \mathcal{L}^{(1,0,\gamma_2)}, \quad \text{for } \gamma_2 \in X_1, \qquad (36a)
$$

$$
\vec{\Gamma}(\gamma_1) = \Gamma^{(2,\gamma_1,0)} \equiv \mathcal{L}^{(2,\gamma_1,0)}, \quad \text{for } \gamma_1 \in X_2. \qquad (36b)
$$

where $\mathcal{L}^{\omega}$ 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\gamma$ denote the usual Lebesgue measure on $X_1 = \Gamma_2$
and $X_2 = \Gamma_1$.

**Proposition 3.5.** Let $\mathbb{P}$ be the probability measure on $X$ with uniform density $\left( |\Gamma_1| + |\Gamma_2| \right)^{-1} d\gamma$.

1. $\mathbb{P}$ is invariant by the groupoid action.

2. The dynamical system $(X, \Gamma(X), \mathbf{t}, \mathbb{P})$ is uniquely ergodic if and only if the lattices $\mathcal{R}_1, \dots, \mathcal{R}_p$ are incommensurate. In this case, we have the Birkhoff property: for any $f \in C(X)$ and $\omega \in X$,

$$
\lim_{r \to \infty} \frac{1}{\#(B_r \cap \mathcal{L}^\omega)} \sum_{\mathbf{a} \in B_r \cap \mathcal{L}^\omega} f(\mathbf{t}_{-\mathbf{a}}\omega) = \int_X f d\mathbb{P}, \quad (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 $\mathfrak{H}_\omega = \ell^2(\Xi^\omega)$ can be seen as:

- elements of $\mathfrak{H}_{\gamma_2} = \ell^2(\mathcal{R}_1) \oplus \ell^2(\gamma_2 + \mathcal{R}_2)$ for $\omega = \gamma_2 \in X_1$,

- or elements of $\mathfrak{H}_{\gamma_1} = \ell^2(\gamma_1 + \mathcal{R}_1) \oplus \ell^2(\mathcal{R}_2)$ for $\omega = \gamma_1 \in X_2$.
Then, given the decomposition (31) of $\Gamma(X)$, it makes sense to write a block decomposition of functions $f \in C^*(\Gamma(X))$ as:

$$f = \begin{bmatrix} f_{11} & f_{12} \\ f_{21} & f_{22} \end{bmatrix} \qquad (38)$$

where

$$ \begin{align} f_{11} &: \quad \vec{\Gamma}_{11} = \Gamma_2 \times \mathcal{R}_1 \to \mathbb{C}, \notag \\ f_{12} &: \quad \vec{\Gamma}_{12} = \mathbb{R}^d \to \mathbb{C}, \tag{39} \\ f_{21} &: \quad \vec{\Gamma}_{21} = \mathbb{R}^d \to \mathbb{C}, \notag \\ f_{22} &: \quad \vec{\Gamma}_{22} = \Gamma_1 \times \mathcal{R}_2 \to \mathbb{C}. \notag \end{align} $$

Note that the decomposition (38) of the tight-binding hopping parameters into intra- and inter-layer terms is naturally the parameterization used in physics [10]. In particular, inter-layer coefficients are usually represented directly as a continuous function of the relative position of the atoms as in (39), see e.g., [9]. Let us now write the * product defining the algebra:

- for $\gamma_2 \in X_1$ and $\mathbf{m} \in \mathcal{R}_1$,

$$ (f * g)_{11}(\gamma_2, \mathbf{m}) = \sum_{\mathbf{m}' \in \mathcal{R}_1} f_{11}(\gamma_2, \mathbf{m}') g_{11}(\gamma_2 - \overline{\mathbf{m}'}^{\mathcal{R}_2}, \mathbf{m} - \mathbf{m}') + \sum_{\mathbf{q}' \in \gamma_2 + \mathcal{R}_2} f_{12}(\mathbf{q}') g_{21}(\mathbf{m} - \mathbf{q}'); \qquad (40a) $$

- for $\mathbf{q} \in \mathbb{R}^2$,

$$ (f * g)_{12}(\mathbf{q}) = \sum_{\mathbf{m}' \in \mathcal{R}_1} f_{11}(\bar{\mathbf{q}}^{\mathcal{R}_2}, \mathbf{m}') g_{12}(\mathbf{q} - \mathbf{m}') + \sum_{\mathbf{n}' \in \mathcal{R}_2} f_{12}(\mathbf{q} - \mathbf{n}') g_{22}(\overline{\mathbf{n}'} - \overline{\mathbf{q}}^{\mathcal{R}_1}, \mathbf{n}'); \qquad (40b) $$

- for $\mathbf{p} \in \mathbb{R}^2$,

$$ (f * g)_{21}(\mathbf{p}) = \sum_{\mathbf{n}' \in \mathcal{R}_2} f_{22}(\overline{\mathbf{p}}^{\mathcal{R}_1}, \mathbf{n}') g_{21}(\mathbf{p} - \mathbf{n}') + \sum_{\mathbf{m}' \in \mathcal{R}_1} f_{21}(\mathbf{p} - \mathbf{m}') g_{11}(\overline{\mathbf{m}'} - \overline{\mathbf{p}}^{\mathcal{R}_2}, -\mathbf{m}'); \qquad (40c) $$

- for $\gamma_1 \in X_2$ and $\mathbf{n} \in \mathcal{R}_2$,

$$ (f * g)_{22}(\gamma_1, \mathbf{n}) = \sum_{\mathbf{n}' \in \mathcal{R}_2} f_{22}(\gamma_1, \mathbf{n}') g_{22}(\gamma_1 - \overline{\mathbf{n}'}^{\mathcal{R}_1}, \mathbf{n} - \mathbf{n}') + \sum_{\mathbf{p}' \in \gamma_1 + \mathcal{R}_1} f_{21}(\mathbf{p}') g_{12}(\mathbf{n} - \mathbf{p}'). \qquad (40d) $$

We can also write the * operation as follows: for $\gamma_2 \in X_1$, $\gamma_1 \in X_2$ and $\mathbf{m} \in \mathcal{R}_1$, $\mathbf{q} \in \gamma_2$, $\mathbf{p} \in \gamma_1$, $\mathbf{n} \in \mathcal{R}_2$:

$$ \begin{align} (f^*)_{11}(\gamma_2, \mathbf{m}) &= \overline{f_{11}(\gamma_2 - \overline{\mathbf{m}}^{\mathcal{R}_2}, -\mathbf{m})}, & (f^*)_{12}(\mathbf{q}) &= \overline{f_{21}(-\mathbf{q})}, \tag{41} \\ (f^*)_{21}(\mathbf{p}) &= \overline{f_{12}(-\mathbf{p})}, & (f^*)_{22}(\gamma_1, \mathbf{n}) &= \overline{f_{22}(\gamma_1 - \overline{\mathbf{n}}^{\mathcal{R}_1}, -\mathbf{n})}. \notag \end{align} $$

Finally, the representation formula writes as follows:
1. for $\gamma_2 \in X_1$ and $\phi = (\phi_1, \phi_2) \in \ell^2(\mathcal{R}_1) \oplus \ell^2(\gamma_2 + \mathcal{R}_2)$, $\pi_{\gamma_2}(f)\phi$ can be decomposed in $\ell^2(\mathcal{R}_1) \oplus \ell^2(\gamma_2 + \mathcal{R}_2)$ as

$$ \left\{ \begin{array}{l} (\pi_{\gamma_2}(f)\phi)_1(\mathbf{m}) = \sum_{\mathbf{m}' \in \mathcal{R}_1} f_{11}(\gamma_2 - \bar{\mathbf{m}}^{\mathcal{R}_2}, \mathbf{m}' - \mathbf{m})\phi(\mathbf{m}') + \sum_{\mathbf{q}' \in \gamma_2 + \mathcal{R}_2} f_{12}(\mathbf{q}' - \mathbf{m})\phi(\mathbf{q}'), \\ (\pi_{\gamma_2}(f)\phi)_2(\mathbf{q}) = \sum_{\mathbf{m}' \in \mathcal{R}_1} f_{21}(\mathbf{m}' - \mathbf{q})\phi(\mathbf{m}') + \sum_{\mathbf{n}' \in \mathcal{R}_2} f_{22}(-\bar{\mathbf{q}}^{\mathcal{R}_1}, \mathbf{n}')\phi(\mathbf{q} + \mathbf{n}'), \\ \text{for all } \mathbf{m} \in \mathcal{R}_1 \text{ and } \mathbf{q} \in \gamma_2 + \mathcal{R}_2; \end{array} \right. \tag{42} $$

2. for $\gamma_1 \in X_2$ and $\phi = (\phi_1, \phi_2) \in \ell^2(\gamma_1 + \mathcal{R}_1) \oplus \ell^2(\mathcal{R}_2)$, $\pi_{\gamma_1}(f)\phi$ can be decomposed in $\ell^2(\gamma_1 + \mathcal{R}_1) \oplus \ell^2(\mathcal{R}_2)$ as

$$ \left\{ \begin{array}{l} (\pi_{\gamma_1}(f)\phi)_1(\mathbf{p}) = \sum_{\mathbf{m}' \in \mathcal{R}_1} f_{11}(-\bar{\mathbf{p}}^{\mathcal{R}_2}, \mathbf{m}')\phi(\mathbf{p} + \mathbf{m}') + \sum_{\mathbf{n}' \in \mathcal{R}_2} f_{12}(\mathbf{n}' - \mathbf{p})\phi(\mathbf{n}'), \\ (\pi_{\gamma_1}(f)\phi)_2(\mathbf{n}) = \sum_{\mathbf{p}' \in \gamma_1 + \mathcal{R}_1} f_{21}(\mathbf{p}' - \mathbf{n})\phi(\mathbf{p}') + \sum_{\mathbf{n}' \in \mathcal{R}_2} f_{22}(\gamma_1 - \bar{\mathbf{n}}^{\mathcal{R}_1}, \mathbf{n}' - \mathbf{n})\phi(\mathbf{n}'), \\ \text{for all } \mathbf{p} \in \gamma_1 + \mathcal{R}_1 \text{ and } \mathbf{n} \in \mathcal{R}_2. \end{array} \right. \tag{43} $$

**Trace per unit volume and integro-differential calculus.** A classical consequence [3] of the ergodicity of $\mathbb{P}$ on $X$ under the action of $\Gamma(X)$, Proposition 3.5, is that we can characterize uniquely $\mathbb{P}$ under the condition that it is a trace per unit volume:

**Proposition 3.6.** When $\mathcal{R}_1$ and $\mathcal{R}_2$ are incommensurate, the invariant, ergodic probability measure $\mathbb{P}$ is uniquely defined as a trace per unit volume in the sense that it satisfies: for any $f \in C^*(\Gamma(X))$ and $\omega \in X$,

$$ \mathcal{T}_{\mathbb{P}}(f) = \lim_{r \to \infty} \frac{1}{\#(B_r \cap \mathcal{L}^{\omega})} \mathrm{Tr} (\pi_{\omega}(f)|_{B_r}). \tag{44} $$

Moreover, for $j = 1, \dots, d$, the derivation defined by (27) extends easily: for $\gamma_2 \in X_1$ and $\mathbf{n} \in \mathcal{R}_1$, $\mathbf{q} \in \gamma_2$, $\mathbf{p} \in \gamma_1$, $\mathbf{m} \in \mathcal{R}_2$,

$$ \begin{aligned} (\partial_j f)_{11}(\gamma_2, \mathbf{n}) &= i n_j f_{11}(\gamma_2, \mathbf{n}), & (\partial_j f)_{12}(\mathbf{q}) &= i q_j f_{12}(\mathbf{q}), \\ (\partial_j f)_{21}(\mathbf{p}) &= i p_j f_{21}(\mathbf{p}), & (\partial_j f)_{22}(\gamma_1, \mathbf{m}) &= i m_j f_{22}(\gamma_1, \mathbf{m}). \end{aligned} \tag{45a} $$

There exists also a derivation for $j = d + 1$

$$ \begin{aligned} (\partial_{d+1} f)_{11}(\gamma_2, \mathbf{n}) &= 0, & (\partial_{d+1} f)_{12}(\mathbf{q}) &= +i f_{12}(\mathbf{q}), \\ (\partial_{d+1} f)_{21}(\mathbf{p}) &= -i f_{21}(\mathbf{p}), & (\partial_{d+1} f)_{22}(\gamma_1, \mathbf{m}) &= 0. \end{aligned} \tag{45b} $$

Note that the derivation operator $\partial_{d+1}$ is bounded on $\mathcal{A}$. Hence, an element $f$ belongs to $\mathcal{C}^N$ if and only if $\|\partial_{\alpha} f\| < \infty$ for all multi-indexes $\alpha$ such that $\sum_{j=1}^d \alpha_j \le N$.

**Remark 3.7.** The extension to multilayer systems of the formalism introduced in the previous sections is straightforward, though somewhat cumbersome.
### 3.3 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 [2], 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 effective 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 $\rho$ as an initial state (i.e. a positive element of the $C^*$-algebra associated with a quantum system), the effective, collision-averaged one-particle time evolution is given by the Liouville equation [21]

$$ \frac{d\rho}{dt} + \mathcal{L}_{h-q\mathcal{E}(t)\cdot\vec{x}}(\rho) = -\frac{\mathbf{1} - \hat{\kappa}^*}{\tau}(\rho), \qquad (46) $$

where

- $h$ is the effective single-electron tight-binding Hamiltonian,

- $\mathcal{E}$ is a constant or time-harmonic spatially uniform electric field with frequency $\hat{\omega}$,

- $\vec{x}$ is the position operator and $\mathbf{1}$ is the identity operator,

- $q$ is the charge of the carrier (hole or electron),

- $\mathcal{L}_{h-q\mathcal{E}(t)\cdot\vec{x}} = i/\hbar [h, \cdot] - \frac{q}{\hbar}\mathcal{E}(t)\cdot\nabla$ is the Liouvillian operator governing the one-particle time evolution in the absence of collisions,

- $\hat{\kappa}^*$ is a scattering-event-averaged collision efficiency operator,

- $\tau$ is the mean collision time of a Poisson process with law $e^{-t/\tau}dt/\tau$ 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 [21]:

$$ \sigma_{ij}(\hat{\omega}) = \frac{q^2}{\hbar^2} \mathcal{T}_{\mathbb{P}} \left( \partial_i h \left[ (\mathbf{1} - \hat{\kappa}^*) / \tau + \mathcal{L}_h - i\hat{\omega} \right]^{-1} \partial_j f_{\beta,\mu}(h) \right), \qquad (47) $$

where $f_{\beta,\mu}(h) = \frac{1}{1+e^{\beta(h-\mu)}}$ is the Fermi-Dirac one-particle density matrix in the grand-canonical equilibrium with chemical potential (Fermi level) $\mu$ and at temperature $T$ such that $\beta = 1/k_B T$.

In the relaxation time approximation (RTA), the effective relaxation operator ($\mathbf{1} - \hat{\kappa}^*)/\tau$ is replaced by $\mathbf{1}/\tau_{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 [2, 3, 19, 21].
# 4 Numerical example

We propose in this section a minimalistic one-dimensional toy model to study incommensurability effects 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 [19], 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 effort is currently ongoing.

## 4.1 Description of the model

Let us consider two parallel one-dimensional crystals with lattice constants $\ell_1$ and $\ell_2$ normalized such that

$$ \ell_1\ell_2 = 1. \qquad (48) $$

Following the notation introduced in Sections 2 and 3, we set the lattices $\mathcal{R}_j = \ell_j\mathbb{Z}$ and the unit cells $\Gamma_j = \mathbb{R}/\mathcal{R}_j = \ell_j\mathbb{T}$ for $j = 1, 2$, where $\mathbb{T} := \mathbb{R}/\mathbb{Z}$ is the one-dimensional 1-periodic torus. The set of all possible configurations is parameterized by the hull $\Omega = \Gamma_1 \times \Gamma_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 \cup X_2$ with $X_1 = \ell_2\mathbb{T}$ and $X_2 = \ell_1\mathbb{T}$.

<figure><img src="image_6.png" /><figcaption>Figure 5: Sketch of the geometry and tight-binding hoppings of our 1D toy model.</figcaption></figure>

We consider the model Hamiltonian generated by an element $h$ of the $C^*$-algebra presented in Section 3.2 and which takes the values

$$ \begin{aligned} h_{11}(\gamma_2, m) &= \begin{cases} 1 & \text{if } |m| = \ell_1, \\ 0 & \text{otherwise,} \end{cases} & \text{for } \gamma_2 \in \ell_2\mathbb{T}, \quad m \in \ell_1\mathbb{Z}, \\ h_{12}(p) &= We^{-\frac{1}{2}\left(\frac{p}{r_0}\right)^2}, & \text{for } p \in \mathbb{R}, \\ h_{21}(q) &= We^{-\frac{1}{2}\left(\frac{q}{r_0}\right)^2}, & \text{for } q \in \mathbb{R}, \\ h_{22}(\gamma_1, n) &= \begin{cases} 1 & \text{if } |n| = \ell_2, \\ 0 & \text{otherwise,} \end{cases} & \text{for } \gamma_1 \in \ell_1\mathbb{T}, \quad n \in \ell_2\mathbb{Z}. \end{aligned} \qquad (49) $$

The corresponding Hamiltonian of the model $H_\omega = \pi_\omega(h)$, 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
depending on the distance between the lattice sites. The two parameters of the model are the maximum amplitude W and the characteristic length r₀ of the inter-chain hopping interaction terms.

## 4.2 Numerical approach

Using this toy model, we aim to showcase the expected effects of continuously varying the lattice constants ratio $\alpha = \ell_2/\ell_1$ in a range around the periodic case of matched chains, $\ell_1 = \ell_2$ i.e. $\alpha = 1$.

### 4.2.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 [12, 19], 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

$$ \ell_1 = \sqrt{\frac{q}{p}}, \quad \ell_2 = \sqrt{\frac{p}{q}} \quad \text{so that} \quad \ell_1 \ell_2 = 1 \quad \text{and} \quad \alpha = \frac{\ell_2}{\ell_1} = \frac{p}{q}. \quad (50) $$

Note that the periodic supercell length is $\sqrt{pq} = p\ell_1 = q\ell_2$. We then scan the range of ratios $1/6 \le \alpha \le 6$ by varying p from $p_{min} = \lfloor \frac{N}{7} \rfloor$ to $p_{max} = \lfloor \frac{6N}{7} \rfloor$. The corresponding tight-binding hamiltonian matrix $H_0^\alpha$ with periodic boundary conditions is then assembled from the $C^*$-algebra element $h^\alpha$ (49) using a cut-off distance of $6r_0$ for the Gaussian inter-layer term, and choosing the configuration $\gamma_1 = \gamma_2 = 0$.

Our two quantities of interest are as follows.

- the density of states, i.e., the spectral measure $d\mu^\alpha(E)$ on $\mathbb{R}$ with support on $\sigma(h)$ is defined by

$$ \int_{\mathbb{R}} \phi(E) d\mu^{\alpha}(E) = \mathcal{T}_{\mathbb{P}}(\phi(h^{\alpha})) \approx \frac{1}{N} \mathrm{Tr}(\phi(H_{0}^{\alpha})), \quad (51) $$

where the test function $\phi$ is analytic in an open neighborhood of $\sigma(h^\alpha)$;

- the conductivity at zero frequency given by the Kubo formula (47), which can be computed more efficiently by introducing the current-current correlation measure $d\mathcal{M}^\alpha$:

$$ \sigma = \left(\frac{e}{\hbar}\right)^2 \iint_{\mathbb{R}^2} \frac{f_{\beta,\mu}(E') - f_{\beta,\mu}(E)}{E - E'} \frac{d\mathcal{M}^{\alpha}(E, E')}{1/\tau_{\mathrm{rel}} - i/\hbar(E - E') - i\hat{\omega}}. \quad (52) $$

The spectral measure $d\mathcal{M}^\alpha$ on $\mathbb{R}^2$ is defined by [20, 21]:

$$ \iint_{\mathbb{R}^2} \phi_1(E)\phi_2(E')d\mathcal{M}^\alpha(E, E') = \mathcal{T}_{\mathbb{P}}(\phi_1(h) \cdot \partial_1 h \cdot \phi_2(h) \cdot \partial_1 h) \quad (53) \\ \approx \frac{1}{N}\mathrm{Tr}(\phi_1(H_0^\alpha) \cdot \widetilde{\partial_1 H_0^\alpha} \cdot \phi_2(H_0^\alpha) \cdot \widetilde{\partial_1 H_0^\alpha}), $$

where the test functions $\phi_1, \phi_2$ are analytic in an open neighborhood of $\sigma(h^\alpha)$. Since $x \mapsto x$ is not $\sqrt{pq}$-periodic, we have used in (53) the approximate periodic differential calculus introduced by Prodan [19], such that the approximate derivation $\widetilde{\partial_1 H_0^\alpha}$ is obtained by:

$$ \widetilde{\partial_1 f}(\omega, x) = i \tilde{\mathcal{X}}(x) f(\omega, x), \quad (54) $$

where $x \mapsto \tilde{\mathcal{X}}(x)$ is $\sqrt{pq}$-periodic and approximates the identity $x \mapsto x$ near zero.
**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 $\mu$ or Boltzmann factor $\beta = 1/k_B T$.

### 4.2.2 Kernel polynomial method

The second ingredient in our calculations is a Chebyshev polynomial expansion, commonly called the Kernel Polynomial Method [26]. The idea is to rescale the Hamiltonian, $H_0^\alpha \to \hat{H}_0^\alpha = (H_0^\alpha - b)/a$ with $a, b$ well chosen so that the spectrum $\sigma(\hat{H}_0^\alpha)$ lies in the energy range $(-1, 1)$.

**Density of States.** The moments of the spectral measure $d\mu^\alpha$ can be computed on the basis of Chebyshev polynomials of the first kind $\{T_m\}_{m \ge 0}$, which form an orthonormal basis of $L^2([-1, 1])$ with respect to the weight function $w(x) = 1/(\pi\sqrt{1-x^2})$. These polynomials obey the recursion relation,

$$ \begin{aligned} T_0(x) &= 1, & T_1(x) &= x, \\ T_{m+1}(x) &= 2xT_m(x) - T_{m-1}(x) & \text{for } m \ge 1. \end{aligned} \qquad (55) $$

From (51), we thus obtain that

$$ \mu_m^\alpha = \int_{\mathbb{R}} T_m \left( \frac{E-b}{a} \right) d\mu^\alpha(E) \approx \frac{1}{N} \mathrm{Tr} \left( T_m(\hat{H}_0^\alpha) \right) = \frac{1}{N} \sum_{j=1}^N T_m(\hat{\lambda}_j^\alpha), \qquad (56) $$

where $\{\hat{\lambda}_j^\alpha\}_{1 \le j \le N}$ is the vector of eigenvalues of $\hat{H}_0^\alpha$. Utilizing the recursion (55), the moments $(\mu_m^\alpha)_{0 \le m \le M}$ can be computed efficiently up to some fixed polynomial degree $M$. Now, assuming that $h^\alpha$ has an absolutely continuous spectrum, the density of states is continuous with respect to the Lebesgue measure, $d\mu^\alpha(E) = \nu^\alpha(E)dE$, we can reconstruct accurately the spectral density $\nu^\alpha$ from the Chebyshev moments [26]:

$$ \nu^\alpha(E) = \frac{1}{\pi\sqrt{a^2 - (E-b)^2}} \left( \mu_0 + 2 \sum_{m=1}^{M} \mu_m^\alpha g_m^M T_m \left( \frac{E-b}{a} \right) \right), \qquad (57) $$

where $g_m^M = [(M-m+1)\cos(\frac{\pi m}{M+1})+\sin(\frac{\pi m}{M+1})\cot(\frac{\pi}{M+1})]/(M+1)$ are the Jackson damping coefficients designed to avoid spurious Gibbs oscillations.

Finally, we note that the values of $\nu^\alpha$ at the particular set of points

$$ x_k = a \cos \left( \frac{\pi(k + 1/2)}{M} \right) + b \quad \text{with } k = 0, \dots, M-1, \qquad (58) $$

coinciding with the abscissas of the Chebyshev-Gauss numerical integration points, can be obtained through a fast cosine transform in $O(M \log M)$ operations since

$$ \gamma_k^\alpha = \pi \sqrt{a^2 - (x_k - b)^2} \nu^\alpha(x_k) = \sum_{m=0}^{M} (1 + \delta_{m,0}) \mu_m^\alpha g_m^M \cos \left( \frac{\pi m (k + 1/2)}{M} \right). $$

This further diminishes the numerical cost of evaluating (57).
**Conductivity.** Let us define the two-dimensional Chebyshev moments of the current-current correlation measure $d\mathcal{M}^\alpha$ following (53):

$$ \mathcal{M}_{mn}^{\alpha} = \iint_{\mathbb{R}^2} T_m\left(\frac{E-b}{a}\right) T_n\left(\frac{E'-b}{a}\right) d\mathcal{M}^{\alpha}(E, E') \\ \approx \frac{1}{N} \text{Tr} \left( T_m(\hat{H}_0^{\alpha}) \cdot \tilde{\partial}_1 \hat{H}_0^{\alpha} \cdot T_n(\hat{H}_0^{\alpha}) \cdot \tilde{\partial}_1 \hat{H}_0^{\alpha} \right). $$

To simplify the computation, let us introduce the diagonal matrix of eigenvalues $\hat{D}^\alpha = \text{Diag}(\lambda_1^\alpha, \dots, \lambda_N^\alpha)$, a unitary matrix of eigenvectors $V^\alpha$ of $\hat{H}_0^\alpha$, and the Hermitian matrix $J^\alpha$ such that

$$ J^{\alpha} = (V^{\alpha})^{*} \tilde{\partial}_{1} \hat{H}_{0}^{\alpha} V^{\alpha} \quad \text{and} \quad \hat{H}_{0}^{\alpha} = V^{\alpha} \hat{D}^{\alpha} (V^{\alpha})^{*}. $$

Then,

$$ \mathcal{M}_{mn}^{\alpha} \approx \frac{1}{N} \text{Tr} \left( T_m(\hat{D}^{\alpha}) \cdot J^{\alpha} \cdot T_n(\hat{D}^{\alpha}) \cdot J^{\alpha} \right) = \frac{1}{N} \sum_{i,j=1}^{N} T_m(\hat{\lambda}_i^{\alpha}) |J_{ij}^{\alpha}|^2 T_n(\hat{\lambda}_j^{\alpha}). \quad (59) $$

The moments $\mathcal{M}_{mn}^\alpha$ can thus be efficiently 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 $\Phi_{\beta,\mu,\tau_{\text{rel}},\hat{\omega}}(E, E')$ the integrand in the right-hand side of (52), we introduce approximate Chebyshev moments $\Phi_{\beta,\mu,\tau_{\text{rel}},\hat{\omega}}^{mn}$ computed with Chebyshev-Gauss integration on the points $x_k$ defined by (58):

$$ \Phi_{\beta, \mu, \tau_{\text{rel}}, \hat{\omega}}^{mn} = \frac{(1 + \delta_{m,0})(1 + \delta_{n,0})}{M^2} \sum_{k,l=0}^{M-1} \Phi_{\beta, \mu, \tau_{\text{rel}}, \hat{\omega}}(b + ax_k, b + ax_l) T_m(x_k) T_n(x_l) \\ \approx (1 + \delta_{m,0})(1 + \delta_{n,0}) \iint_{[-1,1]^2} \frac{\Phi_{\beta, \mu, \tau_{\text{rel}}, \hat{\omega}}(b + ax, b + ay) T_m(x) T_n(y)}{\pi^2 \sqrt{1 - x^2} \sqrt{1 - y^2}} dx dy, $$

so that the function $\Phi_{\beta,\mu,\tau_{\text{rel}}}$ is well approximated by the expansion

$$ \Phi_{\beta, \mu, \tau_{\text{rel}}, \hat{\omega}}(E, E') \approx \sum_{m,n=0}^{M} \Phi_{\beta, \mu, \tau_{\text{rel}}, \hat{\omega}}^{mn} g_m^M g_n^M T_m\left(\frac{E-b}{a}\right) T_n\left(\frac{E'-b}{a}\right). $$

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

$$ \sigma = \iint_{\mathbb{R}^2} \Phi_{\beta, \mu, \tau_{\text{rel}}, \hat{\omega}}(E, E') d\mathcal{M}^{\alpha}(E, E') \approx \sum_{m,n=0}^{M} \Phi_{\beta, \mu, \tau_{\text{rel}}, \hat{\omega}}^{mn} g_m^M g_n^M \mathcal{M}_{mn}^{\alpha} \\ \approx \sum_{m,n=0}^{M} \frac{(1 + \delta_{m,0})(1 + \delta_{n,0})}{M^2} \sum_{k,l=0}^{M-1} \Phi_{\beta, \mu, \tau_{\text{rel}}, \hat{\omega}}(b + ax_k, b + ax_l) T_m(x_k) T_n(x_l) g_m^M g_n^M \mathcal{M}_{mn}^{\alpha}. $$

Finally, we use the relations $T_m(x_k) = \cos(\frac{\pi m(k+1/2)}{M})$ and $T_n(x_l) = \cos(\frac{\pi n(l+1/2)}{M})$ and we exchange the order of summations. The numerical approximation to the conductivity is then given by the quadrature formula

$$ \sigma \approx \frac{1}{M^2} \sum_{k,l=0}^{M-1} \Gamma_{kl}^{\alpha} \Phi_{\beta, \mu, \tau_{\text{rel}}, \hat{\omega}}(ax_k + b, ax_l + b), \quad (60) $$
where the $\Gamma_{kl}^{\alpha}$ are the raw output of the 2D fast cosine transform

$$ \Gamma_{kl}^{\alpha} = \sum_{m,n=0}^{M} (1 + \delta_{m,0})(1 + \delta_{n,0}) \mathcal{M}_{mn}^{\alpha} g_{m}^{M} g_{n}^{M} \cos \left( \frac{\pi m(k + 1/2)}{M} \right) \cos \left( \frac{\pi n(l + 1/2)}{M} \right). $$

In conclusion, formulae (57) and (60) show that the Chebyshev expansion leads to an efficient representation of the spectral measures, relying on fast discrete cosine transforms to compute weights $\gamma_k^\alpha$ and $\Gamma_{kl}^\alpha$ that depend only on the matrix $H_0^\alpha$. The knowledge of these coefficients 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 $\mu, \beta, \tau_{rel}$ or $\hat{\omega}$ for the conductivity.

**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 effective 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 [26].

## 4.3 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 [4]. We plot first the density of states as a function of lattice constants ratio $\alpha$ and 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 [12], 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 $\alpha = \ell_2/\ell_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 ratio $\alpha$ and Fermi level $\mu$ in Figure 7. The same fractal pattern emerges, however a striking difference 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 $\alpha = 1$.

This fan structure is repeated at a number of values of $\alpha$ corresponding 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 $\alpha$.

<table><caption>Table 1: Choice of numerical parameters</caption><thead><tr><th>N</th><th>p<sub>min</sub></th><th>p<sub>max</sub></th><th>M</th><th>W</th><th>r<sub>0</sub></th><th>β</th><th>τ<sub>rel</sub></th><th>ω̂</th></tr></thead><tbody><tr><td>4181</td><td>597</td><td>3583</td><td>1000</td><td>.5</td><td>.25</td><td>250</td><td>250</td><td>0</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 6: Density of states in color scale as a function of energy and lattice constants ratios.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Conductivity in color scale as a function of Fermi level and lattice constants ratios, with parameter choice $W = .5$, $r_0 = .25$, $\beta = \tau = 250$.</figcaption></figure>

Finally, we plot in Figure 8 the conductivity again, which we plot this time as a
function of the integrated density of states

$$n^{\alpha}(E) = \int_{-\infty}^{E} d\mu^{\alpha},$$

where $d\mu^\alpha$ is the density of states measure defined by (51). The function $E \mapsto n^\alpha(E)$ 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 [22], 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 [5], leading to exciting perspectives for the application of the framework presented in this paper to more realistic 2D multilayer systems.

# 5 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.

# References

[1] W. Arveson. A short course on spectral theory, volume 209. Springer Science & Business Media, 2006.

<figure><img src="image_9.png" /><figcaption>Figure 8: Conductivity in color scale as a function of the integrated density of states.</figcaption></figure>
[2] J. Bellissard. Coherent and dissipative transport in aperiodic solids. In *Lecture Notes in Physics*, volume 597, pages 413–486. Springer, 2003.

[3] J. Bellissard, A. van Elst, and H. Schulz-Baldes. The noncommutative geometry of the quantum Hall effect. *Journal of Mathematical Physics*, 35(10):5373–5451, 1994.

[4] J. Bezanson, S. Karpinski, V. B. Shah, and A. Edelman. Julia: A fast dynamic language for technical computing. *arXiv e-print arXiv:1209.5145*, 2012.

[5] Y. Cao, J. Y. Luo, V. Fatemi, S. Fang, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero. Superlattice-induced insulating states and valley-protected orbits in twisted bilayer graphene. *Physical Review Letters*, 117(11):116804, 2016.

[6] S. Carr, D. Massatt, P. Cazeaux, S. Fang, M. Luskin, and E. Kaxiras. Twistronics: Manipulating the electronic properties of two-dimensional layered structures through their twist angle. *arXiv preprint arXiv:1611.00649*, 2016.

[7] P. Cazeaux, M. Luskin, and E. B. Tadmor. Analysis of rippling in incommensurate one-dimensional coupled chains. to appear in *Multiscale Modeling and Simulation*, June 2016.

[8] M. Einsiedler and T. Ward. *Ergodic Theory: with a view towards Number Theory*. Graduate Texts in Mathematics. Springer London, 2010.

[9] S. Fang and E. Kaxiras. Electronic structure theory of weakly interacting bilayers. *Physical Review B*, 93(23):235153, June 2016.

[10] S. Fang, R. Kuate Defo, S. N. Shirodkar, S. Lieu, G. A. Tritsaris, and E. Kaxiras. Ab initio tight-binding hamiltonian for transition metal dichalcogenides. *Physical Review B*, 92(20):205108, November 2015.

[11] A. K. Geim and I. V. Grigorieva. Van der Waals heterostructures. *Nature*, 499(7459):419–25, July 2013.

[12] D. R. Hofstadter. Energy levels and wave functions of bloch electrons in rational and irrational magnetic fields. *Physical Review B*, 14:2239–2249, September 1976.

[13] E. Kaxiras. *Atomic and Electronic Structure of Solids*. Cambridge University Press, Cambridge, 2003.

[14] W. Kohn. Analytic properties of Bloch waves and Wannier functions. *Physical Review*, 115(4):809, 1959.

[15] N. Marzari and D. Vanderbilt. Maximally localized generalized Wannier functions for composite energy bands. *Physical review B*, 56(20):12847, 1997.

[16] D. Massatt, M. Luskin, and C. Ortner. Electronic density of states for incommensurate layers. *arXiv preprint arXiv:1608.01968*, August 2016.

[17] Given $\mathbf{a} \in \mathbb{R}^d$, the translation $T_{\mathbf{a}}$ acts on the space of continuous functions with compact support $C_c(\mathbb{R}^{d+1})$ through $T_{\mathbf{a}}f(\mathbf{x}) = f(T_{-\mathbf{a}}\mathbf{x})$. Therefore it acts on the space $\mathfrak{M}(\mathbb{R}^{d+1})$ of Radon measures through $T_{\mathbf{a}}\mu(f) = \mu(T_{-\mathbf{a}}f)$ whenever $f \in C_c(\mathbb{R}^{d+1})$ and $\mu \in \mathfrak{M}(\mathbb{R}^{d+1})$.
[18] Recall that a cocompact lattice $\mathcal{R}$ of $\mathbb{R}^d$ is a discrete periodic lattice of $\mathbb{R}^d$ with $d$ free vectors (in such a way that $\mathbb{R}^d/\mathcal{R}$ is compact), and that the dual lattice of a cocompact lattice $\mathcal{R}$ is defined by

$$ \mathcal{R}^* = \{\mathbf{k} \in \mathbb{R}^d \mid \mathbf{k} \cdot \mathbf{n} \in \mathbb{Z}, \ \forall \mathbf{n} \in \mathcal{R}\}. $$

When $\mathcal{R}$ is generated by the columns of a matrix $\mathbf{E}$, its dual $\mathcal{R}^*$ is generated by the columns of $(\mathbf{E}^T)^{-1}$, the inverse transpose of $\mathbf{E}$.

[19] E. Prodan. Quantum transport in disordered systems under magnetic fields: A study based on operator algebras. *Applied Mathematics Research eXpress*, 2013(2):176-265, 2013.

[20] E. Prodan and J. Bellissard. Mapping the current-current correlation function near a quantum critical point. *Annals of Physics*, 368:1-15, 2016.

[21] H Schulz-Baldes and J Bellissard. A kinetic theory for quantum transport in aperiodic media. *Journal of Statistical Physics*, 91(5-6):991-1026, 1998.

[22] K. Sugawara, K. Kanetani, T. Sato, and T. Takahashi. Fabrication of Li-intercalated bilayer graphene. *AIP Advances*, 1(2), 2011.

[23] H. Terrones and M. Terrones. Bilayers of transition metal dichalcogenides: Different stackings and heterostructures. *Journal of Materials Research*, 29:373-382, 2014.

[24] G. A. Tritsaris, S. N. Shirodkar, E. Kaxiras, P. Cazeaux, M. Luskin, P. Plecháč, and E. Cancès. Perturbation theory for weakly coupled two-dimensional layers. *Journal of Materials Research*, 31:959-966, 4 2016.

[25] G. H Wannier. The structure of electronic excitation levels in insulating crystals. *Physical Review*, 52(3):191, 1937.

[26] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske. The kernel polynomial method. *Reviews of Modern Physics*, 78:275-306, Mar 2006.