\documentclass[prl,twocolumn,showpacs,superscriptaddress,floatfix,longbibliography]{revtex4-1}
\usepackage{graphicx,amsfonts,amssymb,amsmath,xspace,dsfont}
\usepackage[colorlinks=true,citecolor=green,linkcolor=red,urlcolor=blue]{hyperref}
\usepackage{subfigure}


\usepackage{color}
\definecolor{g}{rgb}{.1,0.4,.1} % {.0,0.7,.5}
\definecolor{b}{rgb}{0,0.2,1}
\definecolor{rouge}{rgb}{0.82,0.,0.}
\definecolor{vert}{rgb}{0.,0.82,0.}
\definecolor{orange}{rgb}{1,0.5,0.}
\definecolor{bleu}{rgb}{0.,0.,0.82}
\definecolor{m}{rgb}{0.82,0.,0.82}
\definecolor{vert2}{rgb}{0.,0.5,0.}
\definecolor{rougeclair}{rgb}{1.0,0.7,0.7}

\newcommand{\red}{\color{rouge}}
\newcommand{\blue}{\color{bleu}}
\newcommand{\green}{\color{vert}}
\newcommand{\orange}{\color{orange}}
\newcommand{\magenta}{\color{m}}
\newcommand{\black}{\color{black}}
\newcommand{\citecolor}{\color{vert2}}

\newcommand{\maH}{\mathcal{H}}
\newcommand{\beq}{\begin{equation}}
\newcommand{\be}{\begin{equation}}
\newcommand{\beqn}{\begin{eqnarray}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\eeqn}{\end{eqnarray}}
\newcommand{\nn}{\nonumber}
\newcommand{\tr}{\textrm}
\newcommand{\q}{{\bf q}}
\newcommand{\ep}{{\epsilon}}
\renewcommand{\k}{{\bf k}}
\renewcommand{\r}{{\bf r}}
\newcommand{\K}{{\bf K}}
\renewcommand{\v}{{\bf v}}
\newcommand{\D}{{\bf D}}
\newcommand{\bQ}{{\bf Q}}
\renewcommand{\a}{{\bf a}}
\newcommand{\G}{{\bf G}}
\newcommand{\bem}{\begin{pmatrix}}
\newcommand{\eem}{\end{pmatrix}}
\newcommand{\hb}{\hbar}
\newcommand{\8}{\infty}
\newcommand{\f}{\frac}
\newcommand{\jn}[1]{\textcolor{red}{#1}}

\newlength{\ldag}
\settowidth{\ldag}{\scriptsize $\dagger$}
\newcommand{\cra}{c^\dagger}
\newcommand{\ana}{{c^{\phantom\dagger}\hspace{-\ldag}}}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
\begin{document}

\title{Hofstadter butterfly of a quasicrystal}


\author{Jean-No\"el Fuchs}
\email{fuchs@lptmc.jussieu.fr}
\affiliation{Laboratoire de Physique Th\'eorique de la Mati\`ere Condens\'ee, CNRS UMR 7600, Universit\'e Pierre et Marie Curie, 4 Place Jussieu, 75252 Paris Cedex 05, France}
\affiliation{Laboratoire de Physique des Solides, CNRS UMR 8502, Universit\'e Paris-Sud, 91405 Orsay, France}

\author{Julien Vidal}
\email{vidal@lptmc.jussieu.fr}
\affiliation{Laboratoire de Physique Th\'eorique de la Mati\`ere Condens\'ee, CNRS UMR 7600, Universit\'e Pierre et Marie Curie, 4 Place Jussieu, 75252 Paris Cedex 05, France}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}

The energy spectrum of a tight-binding Hamiltonian is studied for the two-dimensional quasiperiodic Rauzy tiling in a perpendicular magnetic field. This spectrum known as Hofstadter butterfly displays a very rich pattern of bulk gaps that can be labeled by four integers, instead of two integers for periodic systems. The role of phason-flip disorder is also investigated in order to extract the genuinely quasiperiodic property. This geometric disorder is found to only preserve main quantum Hall gaps.

\end{abstract}

\pacs{??}

\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\textit{Introduction.--- }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
{\bf Quasicrystals are nonperiodic solids that nevertheless feature long-range order \cite{Shechtman84}. In reciprocal space, this order is characterized by resolution-limited Bragg peaks in the diffraction pattern. In real space, the order is related to the nonperiodic repetitivity of local environments. 
After the initial burst of interest in the 80's and 90's, there is a recent revival of activities on the transport properties of quasicrystals mainly due to new experimental realizations of artificial quasiperiodic lattices with phonons, cold atoms \cite{Verkerk99}, photons \cite{Krauss12}, polaritons \cite{Tanese14} or microwaves \cite{Vignolo16}. These new tools allows one to experimentally address questions that were impossible to settle in solid state (metallic alloys) quasicrystals, such as the labeling of gaps and the measure of the energy spectrum or the nature of the eigenstates wave functions. 

A current trend consists in connecting the properties of quasicrystals to that of topological insulators (see, e.g., Refs.~\cite{Krauss12,Tran15}). In one dimension, integers used to label the gaps of the Fibonacci chain have been related to Chern numbers \cite{Levy15}. Although, there still is a controversy on whether one-dimensional quasicrystals are truly topological insulators themselves or whether this is a property of a family of such quasicrystals \cite{Madsen13}, just as quantized adiabatic pumping in one dimension is actually related to the integer quantum Hall effect (IQHE) in two dimensions. 

In the present theoretical work, we use the orbital response of electrons to a magnetic field to probe the quasiperiodic order in a simple two-dimensional tight-binding model. In particular we study the energy spectrum in a magnetic field, i.e., the Hofstadter butterfly \cite{Hofstadter76}. We show that the gap labeling in this case involves both the topological numbers of the IQHE and integers related to the irrationals used in the construction of quasicrystals.
}
\medskip

%
%%%%%%%%%%%
\begin{figure}[t]
\includegraphics[width=\columnwidth]{./figures/flip/flip.pdf}
\caption{A piece of the isometric Rauzy tiling (left) transformed into a random tiling (right) via a random sequence of phason flips (middle). The red dot indicates the center of the inversion symmetry which is broken after flips.}
\label{fig:tiling}
\end{figure}
%%%%%%%%%%%
%
%

%
%
%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure*}[t]
\centering
\includegraphics[width=1\textwidth]{./figures/Butterflies/Butterflies.pdf}
\caption{(Color online) Left: Hofstadter butterflies of the isometric Rauzy tiling on a torus with $R_{15}=5768$ sites (left) and of a random tiling obtained after $2.10^7$ flips (right). Arrows indicate some fluxes around which Landau levels are found with a nonvanishing Hall conductivity.}
\label{fig:pappbc}
\end{figure*} 
%%%%%%%%%%%%%%%%%%%%%%%
%
%

%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\textit{Rauzy and random tilings.---}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
Rauzy tilings can be seen as generalizations of the Fibonacci chain to higher dimensions \cite{Vidal01}. These codimension-1 quasicrystalline tilings are built using the celebrated cut-and-project method \cite{Kalugin85, Duneau85, Elser85}. In the following, we consider the two-dimensional Rauzy tiling and its approximants whose construction  is based on the Tribonacci sequence defined as 
%
%
%%%%%%%%%%%
\be 
R_{n+1}=R_n+R_{n-1}+R_{n-2}, \: \forall n\in \mathbb{N},
\ee 
%%%%%%%%%%%
%
%
with \mbox{$R_0=R_1=1$}, and $R_2=2$. The order-$k$ approximant contains $R_{k+1}$ sites and, after a proper ordering of the sites (according to their position in the perpendicular space in the cut-and-project construction), its connectivity matrix has a Toeplitz-like structure with bands starting at positions ($R_{k-2}, R_{k-1}, R_{k}$) \cite{Vidal01}. This rhombus tiling contains 3, 4, and 5-fold coordinated sites. In the quasiperiodic limit, their densities are given  by \mbox{$\rho_3=\rho_5=2\,\theta^{-3}$} and \mbox{$\rho_4=1-\rho_3-\rho_5$}, where $\theta\simeq 1.839$ is the so-called \mbox{Tribonacci} constant defined as the Pisot root of the equation \mbox{$x^3=x^2+x+1$}. As any codimension-1 tiling, Rauzy tiling approximants only possess an inversion symmetry associated with the center of the one-dimensional acceptance zone. 

In its original construction \cite{Vidal01}, the Rauzy tiling has three different types of tiles (corresponding to the projections of the cubic-lattice faces onto the parallel space) with incommensurate areas. However, one can deform the tiling such that all areas become identical by changing the projection direction in the cut-and-project algorithm.  This isometric version of the Rauzy tiling displayed in Fig.~\ref{fig:tiling} is especially well-suited to the problem under study \cite{Vidal04} (see below).  Moreover, we will also pay attention to a structural disorder induced by phason flips which consist in locally changing neighbors of 3-fold coordinated sites as depicted in Fig.~\ref{fig:tiling}. This does not change the number of links in each direction. As argued in Ref.~\cite{Destainville02} one needs to perform about $N^2/2$ flips to fully disorder a tiling with $N$ sites. After such a rearrangement of links, one obtains a random tiling with 3, 4, 5, and 6-fold coordinated sites.
% with densities that depend on the initial tiling. Starting from the Rauzy tiling, one gets equilibrium densities approximately given by $\rho_3 \simeq 0.19$, $\rho_4\simeq 0.63$, $\rho_5 \simeq 0.17$ and $\rho_6 \simeq 0.01$. 



\medskip


%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\textit{Model and symmetries.---}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
For simplicity, we consider a single-orbital tight-binding Hamiltonian 
%
%
%%%%%%%%%%%%
\be 
H=-\sum_{\langle i,j \rangle} t_{ij} \: |i\rangle \langle j|
%{c^\dagger_i} {{c^{\phantom\dagger}\hspace{-\ldag}_j}}
,
\label{eq:ham}
\ee
%%%%%%%%%%%%
%
%
%where ${c^\dagger_i}$ $({{c^{\phantom\dagger}\hspace{-\ldag}_i}})$ are creation (annihilation) fermionic operator on site $i$ and  
where the sum is performed over nearest-neighbor sites. When a magnetic field $\boldsymbol{B}$ perpendicular to the tiling is introduced, the hopping term from site $i$ to site $j$ is modified according to the Peierls substitution \mbox{$t_{ij} \rightarrow t_{ij} \, {\rm e}^{-\frac{2 {\rm i} \pi}{\phi_0} \int_{i}^j d\boldsymbol{r}\cdot \boldsymbol{A}(\boldsymbol{r})}$} where $\boldsymbol{A}$ is a vector potential such that $\boldsymbol{B}=\boldsymbol{\nabla}\times \boldsymbol{A}$. In the following, we use the following units: $t_{ij}=1$, $\hbar=1$, $e=-2\pi$ so that the flux quantum $\phi_0=h/|e|=1$, and the nearest neighbor distance $a=1$. We also introduce the reduced flux per  plaquette $f=\phi / \phi_0=\pm |\boldsymbol{B}| {\mathcal A}$, where ${\mathcal A}=\frac{\sqrt{3}}{2} a^2$ is the area of an elementary rhombus.

Since ${\mathcal A}$ is the same for all rhombi, the spectrum of $H$ is periodic with $f$ (at least for open boundary conditions). Consequently, one can restrict the study to $f \in [0,1]$. The spectrum is also obviously unchanged when the field direction is reversed  ($f \leftrightarrow -f$). In addition, since the lattice is bipartite, the spectrum is  symmetric with respect to $0$. For periodic boundary conditions, this symmetry is broken due to odd cycles encircling the torus that destroy bipartiteness. Similarly, in the presence of a magnetic field, fluxes are present in the torus and destroy the periodicity with $f$.  However, these two symmetry-breaking effects become negligible in the thermodynamic limit. 

\medskip


%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\textit{Boundary conditions and gauge choice.---}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
Since the pioneering work of Hofstadter on the square lattice \cite{Hofstadter76}, the spectrum of $H$ as a function of $f$, dubbed ``Hofstadter butterfly", has been analyzed for many periodic two-dimensional systems (triangular lattice \cite{Claro79}, honeycomb lattice \cite{Rammal85}, flat-band lattices \cite{Aoki96}, dice lattice \cite{Vidal98}, kagom\'e lattice \cite{Xiao03},...) unveiling very rich features. The simplicity of these structures allows one to study the butterfly directly in the thermodynamical limit using suitable choices for the vector \mbox{potential $\boldsymbol{A}$}. 
For quasiperiodic systems, one needs to consider a finite size system and, for any gauge choice, two problems arise. First, the incommensurability of tile areas breaks the periodicity of the butterfly with $f$, as originally discussed in the Penrose lattice \cite{Hatakeyama89}. Second, as for any other system, if one considers open boundary conditions, edge states prevent one from identifying the bulk gaps as discussed in Refs.~\cite{Vidal04,Tran15}. 

In the present work, we solve these two issues by (i) deforming the tiling to deal with identical tile areas (see discussion above) and (ii) by considering periodic boundary conditions. This latter condition can always be fulfilled at the price of restrictions on the accessible reduced flux $f$ (see supplementary material). Indeed, the total flux through the system must be an integer \cite{Zak64b,LL9}. Thus, for a tiling with $N$ identical plaquettes on a torus, $f$ must be a multiple of $1/N$. Note that, in the isometric Rauzy tiling on a torus, the number of plaquettes equals the number of sites. 
In the following, we consider a unit cell of an approximant with periodic boundary conditions and perform numerical diagonalization of $H$ for systems up to $R_{19} = 66012$ sites.

\medskip

%
%
%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{./figures/Wannier/Wannier2.pdf}
\caption{(Color online) Wannier diagram of the isometric Rauzy tiling on a torus with $R_{16}=10609$ sites. Red lines highlight main IQHE gaps $(\nu,0,0,0)$,  green lines indicate some \mbox{($\nu=0$)-gaps}, and magenta lines illustrate gaps $(\nu,1-\nu,0,0)$. Symmetries allow one to restrict the relevant range of  $f$ and $\mathcal{N}$ to $[0,1/2]$. Only gaps larger than $2.10^{-2}$ are shown. This diagram and the gap labeling are unchanged for larger approximants.}
\label{fig:wannier}
\end{figure} 
%%%%%%%%%%%%%%%%%%%%%%%
%
%

%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\textit{Hofstadter butterfly and Wannier diagram---}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
The zero-field energy spectrum of the Rauzy tiling has been discussed in Refs.~\cite{Jagannathan01,Triozon02}. Contrary to the Penrose tiling \cite{Zijlstra00} or the octagonal tiling \cite{Zijlstra04}, the spectrum of the pure hopping model (\ref{eq:ham}) is gapless (see supp. mat.). In the presence of a magnetic field, the spectrum has been computed for open boundary conditions \cite{Vidal04,Tran15} but,  despite some efforts to get rid of them, edge states fill the bulk gaps emerging for nonvanishing fields so that it is impossible to analyze the nontrivial characteristics butterfly. By contrast, as shown in Fig.~\ref{fig:pappbc} (left), a very rich gap structure is unveiled when considering the system on a torus. We emphasize that all gaps visible are stable when increasing the order of the approximant so that, up to the image resolution, this butterfly should be considered as the one of the (infinite) quasiperiodic isometric Rauzy tiling.  As in most Hofstadter butterflies, one observes the presence of Landau levels arising from band edges separated by IQHE gaps. As usual, these levels are broadened when the system is disordered. Nevertheless, as can be seen in Fig.~\ref{fig:pappbc} (right), the phason-flip disorder is sufficiently weak to preserve these IQHE gaps while destroying the fine structure. 

To go beyond and in the absence of an exact mapping of the Schr\"odinger equation onto a Harper-like \mbox{equation \cite{Hofstadter76}}, we focus on the analysis of gap structure. To this aim, we compute the so-called Wannier diagram \cite{Wannier78} obtained by plotting, for an energy $E$ inside a gap, the normalized integrated  density of states $\mathcal{N}(E,f)$, i.e., the number of levels below $E$ divided by the total number of levels, as a function of $f$ (see Fig.~\ref{fig:wannier}). 
\medskip

%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\textit{Gap labeling.---}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
In the quasiperiodic limit, any gap can be labeled with four integers $(\nu, u, v, w)$. Indeed, integrating
the St$\check{\rm{r}}$eda formula~\cite{Streda82} for the Hall conductivity at energy $E$ inside a gap
%
%
%%%%%%%%%%%%%%%%%%%%%%%
\beq
\sigma =  e  \frac{\partial \mathcal{N}(E,f)}{\partial  \phi}= -\frac{e^2}{h} \nu,
\eeq 
%%%%%%%%%%%%%%%%%%%%%%%
%
%
one finds $\mathcal{N}(E,f)=\nu f + \mathcal{N}(E,0)$, where $\nu$ is an (topologically-invariant) integer \cite{Thouless82,Niu85}. This linear dependence is directly observed in Fig.~\ref{fig:wannier}. For open boundary conditions, $\nu$ counts the number of edge states, as recently discussed in Ref.~\cite{Tran15} for the Rauzy tiling. Since (i) the normalized integrated density of states is a multiple of $1/R_{k+1}$ for the order-$k$ approximant and (ii) three consecutive Tribonacci numbers are coprime integers, B\'ezout's identity guarantees that there exists an integer triplet $(u,v,w)$ such that
%
%
%%%%%%%%%%%%%%%%%%%%%%%
\beq
 \mathcal{N}(E,0)=u\frac{R_{k-1}}{R_{k+1}}+v\frac{R_{k}}{R_{k+1}}+w.
\eeq
%%%%%%%%%%%%%%%%%%%%%%%
%
%
Although we have no rigourous proof, we found that $(u,v,w)$ do not depend on $k$. Thus, in the quasiperiodic limit ($k\rightarrow \infty$), one straightforwardly gets: 
%
%
%%%%%%%%%%%%%%%%%%%%%%%
\beq
\mathcal{N}(E,f)=\nu \, f + u\, \theta^{-2} + v \, \theta^{-1} + w.
\label{eq:conj}
\eeq 
%%%%%%%%%%%%%%%%%%%%%%%
%
Note that this result could certainly be obtained, in a more sophisticated way, using the gap-labeling \mbox{theorem \cite{Baake_book00}}. Moreover, it is clear that this line of reasoning can be applied to any other cut and project tiling built on a Pisot number (Penrose, octagonal, icosahedral,...). 

The Wannier diagram reveals the existence of essentially different types of gaps associated either to quasiperiodic order (destroyed by flips) or to IQHE (preserved by flips). Indeed, as can be seen in Fig.~\ref{fig:pappbc}, only IQHE gaps ($u=v=0$) \cite{Thouless82} survive disorder. Gaps identified in Ref.~\cite{Tran15} belong to this family (see red lines in Fig.~\ref{fig:wannier}). However, in the quasiperiodic case, many other gaps are observed in the diagram (see Fig.~\ref{fig:wannier}). Among them, one may distinguish the one existing at $f=1/2$ where time-reversal symmetry implies a vanishing Hall conductivity. Such gaps are also found  for many other fluxes (see green lines in Fig.~\ref{fig:wannier}), but they might be considered as a curiosity. Indeed, to our knowledge, \mbox{($\nu=0$)-gaps} have only been observed so far in the Lieb \cite{Goldman11} and in the dice lattice \cite{Vidal98}. A general existence condition of these gaps is still lacking. 

Finally, there exists some irrational fluxes in the vicinity of which one finds gaps and Landau levels that are destroyed by disorder but have a nonvanishing Hall conductivity ($\nu\neq 0$). These fluxes are thus deeply related to quasiperiodic order and play a role similar to the rational fluxes in the Hofstadter butterfly of periodic crystals \cite{Hofstadter76}.  As for the gaps, one can label these fluxes in the thermodynamical limit according to $f=p\, \theta^{-2} + q \, \theta^{-1} + r$. These special irrational fluxes are also local minima of the ground-state energy as a function $f$ as for rational fluxes in periodic systems (see Ref.~\cite{Pannetier84} for an experimental observation of this phenomenon in the square lattice). 
For illustration, we show three examples in Fig.~\ref{fig:pappbc} corresponding to $(p,q,r)=(1,0,0), (0,1,0)$, and $(-1,-1,1)$. The Landau-level fan for $f=\theta^{-2}$ is displayed in {Fig.~\ref{fig:wannier} (magenta lines). At this stage, we cannot prove that all possible sets of $(p,q,r)$ give rise to such singularities but we conjecture it is so. 
 
\medskip



%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\textit{Landau levels and effective mass.---} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
To better characterize the IQHE gaps, let us focus on Landau levels that arise from band edges (see Fig.~\ref{fig:pappbc}). In the zero-flux limit, the excitation energy of the $n$-th Landau level $\Delta E_n$ can be well fitted by 
%
%
%%%%%%%%%%%%%
\be
\Delta E_n=\hbar \frac{|e {\bf B}|}{m}\left(n+\frac{1}{2}\right)=\frac{4\pi f}{\sqrt{3} \: m}\left(n+\frac{1}{2}\right), \:  \forall n\in \mathbb{N},
\label{eq:LL}
\ee
%%%%%%%%%%%%%
%
%
where $1/m= 1.957(2)$ is the inverse effective mass of the electron. As expected, for the order-$k$ approximant and a given flux $f$, the degeneracy of each Landau level is given by $R_{k+1} \: f$. However, when $f$ increases, the degeneracy of these levels is lifted since lattice effects lead to a broadening as discussed in Ref.~\cite{Zak64} for crystals. 

There are several ways to understand the surprising emergence of an effective mass in nonperiodic systems (see Supp. Mat.). One possibility is to consider an infinite approximant structure with $R_{k+1}$ sites per unit cell and to compute, for $f=0$, the inverse effective mass tensor $\alpha$ of the lowest-energy band. Practically, one diagonalizes the Bloch Hamiltonian $H(\boldsymbol{k})=e^{-i\boldsymbol{k}\cdot \hat{\boldsymbol{r}}}H e^{i\boldsymbol{k}\cdot \hat{\boldsymbol{r}}}$, where $\boldsymbol{k}$ is the Bloch wave vector and $\hat{\boldsymbol{r}}$ is the position operator. One then expands the dispersion relation of the lowest-energy band in the vicinity of $\boldsymbol{k}=(0,0)$: 
\be
\epsilon (\boldsymbol{k})\simeq \epsilon_0+\frac{1}{2}\alpha_{ij}k_i k_j,
\ee
%
%
%
%%%%%%%%%%%
%\beq
%\underline{\alpha}=\left(\begin{array}{cc}1.61702&0.0803396\\0.0803396&2.37328\end{array} \right).
%\label{alpha}
%\eeq
%%%%%%%%%%%
%
%
where $\epsilon_0=-4.115008(1)$ is the ground-state energy. Denoting $\alpha_1$ and $\alpha_2$ the eigenvalues of $\alpha$, 
the average inverse effective mass is then given by \mbox{$1/{m_\text{T}}=\sqrt{\alpha_1 \alpha_2}=1.95735(1)$}. This effective mass  is related to the Thouless conductance \cite{Edwards72,Abrahams79} for the lowest-energy band through $g \sim 1/{m_\text{T}}$ . 
Actually, it may appear fortuitous that ${m_\text{T}}$ matches $m$ so well as it is computed from the curvature of a single band whereas Landau levels are built from $R_{k+1} \: f$ bands. This result is due to a finite stiffness of the ground-state energy with respect to boundary conditions (Thouless energy \cite{Edwards72}). This finite stiffness stems from the extended nature of the ground state that we have explicitly checked. 

In the disordered case, the situation is different. As already mentioned, Landau levels broaden so that, even in the zero-flux limit, a precise determination of the effective mass is harder. On the one hand, we get $1/m\simeq 2$ by a bruteforce fit of the Landau-level slope. On the other hand, for $f=0$ and since the system is disordered, one expects all states to be localized \cite{Abrahams79} and we indeed find that $1/{m_\text{T}} \propto {\rm e}^{-L/\xi}$  when increasing the linear system size $L$ and with a localization length $\xi \sim 30 \: a $. From that respect, phason disorder should be considered as a weak disorder. When the cyclotron radius is smaller than the localization length, i.e., $f \ll \xi^{-2}$, energy levels are insensitive to the magnetic field. In the opposite case, broad Landau levels show up as can be seen in Fig.~\ref{fig:pappbc}.  

\medskip

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\textit{Conclusion.---}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
{\bf In the present work, we obtained bulk Hofstadter butterflies for the quasiperiodic and the disordered Rauzy tiling. The clean butterfly has a complicated structure involving several types of gaps related either to the magnetic field (IQHE) or to quasiperiodic order or to both. Random phason flips are shown to disorder the tiling and to wash out the inner structure of the butterfly: only main IQHE gaps remain. For the clean Rauzy tiling, we have carefully checked that our results are size converged, which seems to be related to the delocalized nature of eigenstates in the Rauzy tiling.  For future work, it would therefore be interesting to test these ideas higher-codimension quasicrystals such as the octagonal or the Penrose tiling, that have critical eigenstates \cite{Kalugin14}. 

We have found that all gaps in the clean butterfly are labeled by four integers. This is in contrast to periodic systems, for which only two integers are required \cite{Dana85}. The number of integers needed to label gaps depends on the nature of the irrational involved in the construction of the quasicrystal. We also provide a first example of gap labeling involving both IQHE and the quasiperiodicity. At the moment, only integers related to IQHE have a clear topological interpretation \cite{Thouless82,Niu85}.

Many open questions are raised by comparing the butterfly of a quasicrystal and that of a periodic system. For example, why do gaps close at finite fluxes instead of crossing the whole Wannier diagram? Why is an obvious self-similarity nor present in the butterfly? Is its spectral measure finite?}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\acknowledgments

We thank M. Ullmo for her contributions at an early stage of this work, B. Dou\c{c}ot for pointing out the importance of magnetic translations and E. Akkermans, A. Jagannathan, P. Kalugin, J. Kellendonk, N. Mac\'e, G. Montambaux, R. Mosseri, F. Pi\'echon and  A. Soret for stimulating discussions.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibliography{./biblio_Rauzy}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
%
\appendix
%
\section*{Supplementary material}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\section*{Periodic boundary conditions and magnetic field: magnetic translations}
\label{mt}
A difficulty in computing the energy spectrum of a quasiperiodic (or random) tiling in the presence of a perpendicular magnetic field is that one has to work with a finite-size system. Indeed, one cannot use Bloch's theorem, as in the case of a periodic lattice, in order to directly work in the thermodynamic limit. Working with a finite-size system, one has to make a choice for boundary conditions. Open boundary conditions are useful in the sense that any magnetic field is possible. But one drawback is that bulk levels are mixed with edge levels. In the present work, we are interested in bulk properties and, in particular, we want to clearly identify bulk gaps. We therefore want to impose periodic boundary conditions.

A standard approach is to make a gauge choice for the vector potential (such as Landau gauge) and then to try to impose that the Peierls phase matches the periodic boundary conditions. This is actually very inefficient and usually strongly restricts the allowed values of magnetic fluxes that are allowed \cite{Tran15}. In particular, it does not allow one to compute the equivalent of a Hofstadter butterfly.

However, there is a general result known as magnetic translations that can help us \cite{Zak64,LL9}. Magnetic-translation operator are a generalization of the usual translation operators when a magnetic field is present. Indeed, when the hopping amplitudes of a tight-binding Hamiltonian are dressed with Peierls phases, the resulting Hamiltonian $H$ no longer commutes with translation operators $T_{\boldsymbol{a}_j}$ where $\boldsymbol{a}_1$ and $\boldsymbol{a}_2$ are unit cell vectors of an approximant (the total area of the sample being $|\boldsymbol{a}_1\times \boldsymbol{a}_2|$). This is due to the appearance of the vector potential $\boldsymbol{A}$, which is non-uniform even when the magnetic field is. However, $H$ still commutes with magnetic translation operators 
\be
\mathcal{T}_{\boldsymbol{a}_j}=e^{i2\pi \chi_{\boldsymbol{a}_j}(\hat{\boldsymbol{r}})}T_{\boldsymbol{a}_j},
\ee
with 
\be
\chi_{\boldsymbol{a}_j}(\boldsymbol{r})=\int^{\boldsymbol{r}}_{\boldsymbol{0}} d\boldsymbol{r}' \cdot [\boldsymbol{A}(\boldsymbol{r}'-\boldsymbol{a}_j)-\boldsymbol{A}(\boldsymbol{r}')],
\ee
which are product of a gauge transformation $e^{i2\pi \chi_{\boldsymbol{a}_j} (\hat{\boldsymbol{r}})}$ and of a translation operator $T_{\boldsymbol{a}_j}$, where $\hat{\boldsymbol{r}}$ denotes the position operator. Magnetic translation operators along the two directions $\boldsymbol{a}_1$ and $\boldsymbol{a}_2$ do no commute in general. They only commute if the total magnetic flux across the sample $N f$ is a multiple of the flux quantum $\phi_0$. Therefore, $f$ must be a multiple of $1/N$ where $N$ is the number of elementary plaquettes. For simplicity, we assume here that all tiles have the same area.

A concrete implementation of the ``magnetic translation trick'' is to make a specific gauge choice, to compute the Peierls phase on every link of the open boundary system and then to include both a Peierls phase and the gauge transformation phase on the stitches needed to realize periodic boundary conditions, i.e., on the links needed to glue the open system into a torus.

One way of checking the constraint on the flux is to compute the trace of $H^4$, which is related to closed paths of length 4 on the tiling. This gives
\be
\text{Tr }H^4 =  (N-1)\cos[2\pi f] + \cos[(1-N)2\pi f] + \text{const},
\ee
The two first terms mean that there are $N-1$ tiles encircling a flux $f$ and one tile that encircles a flux $(N-1)f$. The last term is independent on $f$ and is due to self retracing paths that do not encircle any flux. One sees that for an arbitrary $f$, all tiles but one are threaded by the same flux. However when $N f$ is an integer $p$ then $(1-N) f = f =\frac{p}{N}$ modulo $1$ and this last tile is equivalent to all others. This explains the restriction to $f=\frac{p}{N}$.

{\bf It is important to realize that this trick can be implemented in any gauge, in contrast to the claim made in \cite{Werner15}. We have explicitly checked that by implementing a family of Landau gauges
\be
\boldsymbol{A}=B(x\cos \theta_g  + y \sin \theta_g)(\cos \theta_g\boldsymbol{u}_y-  \sin \theta_g\boldsymbol{u}_x),
\ee
parametrized by an angle $\theta_g$. A subtle point is that varying $\theta_g$ is not rigorously speaking a pure gauge transformation since, on a finite-size system with periodic boundary conditions, there is not only a magnetic flux in each tile but also on large paths that encircle the torus. Whereas the former flux is independent on $\theta_g$, the latter depends on it. Therefore the energy spectrum is not strictly speaking independent on $\theta_g$. However, this effect is small and vanishes in the thermodynamic limit.
}

\section*{Zero-field  density of states and effective band edge mass}
\label{ebem}

In this section, we present an alternative way of defining an effective band edge mass for a tiling. To this aim, we start by briefly discussing the zero-field thermodynamic density of states 
\be
\rho(\mu,T)=\sum_\alpha \frac{1}{4T}\text{sech}^2\frac{E_\alpha-\mu}{2T},
\ee 
plotted in Fig.~\ref{fig:dos}, where $T$ is the temperature, $\mu$ the chemical potential and $\{E_\alpha\}$ are the energy eigenvalues. Here, we set the Bolzmann constant $k_{\rm B}=1$. Temperature is used to smoothen the density of states and corresponds to a box width of $\Delta E \simeq 3.53 \: T$. 

At high temperatures,  clean and disordered cases coincide for all chemical potentials (see Fig.~\ref{fig:dos}). At low temperatures, the agreement remains good (and the density of states smooth) near band edges, but there is a strong disagreement near the band center. In the disordered case, the density of states is smooth,  apart from a zero-energy delta peak corresponding to less than 1\% of very localized states around 6-fold coordinated sites (see Ref.~\cite{Sutherland86} for a description of these states in the dice lattice). In the clean case, the low-temperature density of states has a lot of structure related to the presence of pseudo-gaps. 
%%
%%
%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\includegraphics[width=\columnwidth]{./figures/DOS/DOS.pdf}
\end{center}
\caption{Thermodynamic density of states per site $\rho(\mu,T)/N$ at temperature $T$ as a function of the chemical potential $\mu$ for  (a) the Rauzy tiling; (b) the random tiling. Results are size-converged.
}
\label{fig:dos}
\end{figure}
%%%%%%%%%%%%
%%
%%
%


In the main text, we have defined an effective (band-edge) mass from the band structure. An alternative way consists in fitting the smoothed band-edge normalized integrated density of states in zero-field $\mathcal{N}(E,f=0)$ assuming a parabolic band edge. This amounts to write 
%
%
%%%%%%%%%%%%%%%%%%%%
\be
\mathcal{N}(E,f=0)=\frac{m_\rho}{2\pi}\frac{\sqrt{3}}{2} (E-\epsilon_0),
\label{eq:1}
\ee
%%%%%%%%%%%%%%%%%%%%
%
% 
where the zero-field ground-state energy  in the clean case is \mbox{$\epsilon_0=-4.115008(1)$} whereas $\epsilon_0=-4.08(1)$ in the disordered case. Fitting the effective mass with this expression, one gets $1/{m_\rho}= 1.95(1)$ in the clean case and $1/{m_\rho}= 1.9(1)$ in the disordered case. However, $m_\rho$ and $m$ discussed in the main text are identified through the Onsager semiclassical quantization of closed cyclotron orbits~\cite{Onsager52}
\be
\mathcal{N}(E_n,f=0)=\left(n+\frac{1}{2}\right)f.
\label{eq:2}
\ee 
%
Using Eqs.~(\ref{eq:1}-\ref{eq:2}), one indeed finds
\be
E_n-\epsilon_0=\frac{4\pi f}{\sqrt{3} \: m_\rho}\left(n+\frac{1}{2}\right),  \: \forall n\in \mathbb{N},
\ee
which is similar to Eq.~(\ref{eq:LL}) provided $m_\rho=m$. We thus have two independent ways to compute this effective mass: a direct fit of the Landau level that gives $m$, and a fit of $\mathcal{N}(E,f=0)$ according to Eq.~(\ref{eq:1}) that gives $m_\rho$. If both approaches are in good agreement for the clean case, results for the disordered case are less precise. A better analysis would require an average over a large number of disorder configurations but this is beyond the scope of the present work. 

A third approach to compute this effective mass is discussed in the main text. It relies on a quadratic expansion of the lowest-energy band near its minimum and gives an effective mass $m_{\rm T}$ related to the Thouless conductance. For the Rauzy tiling, the effective mass tensor $\alpha$ has two different eigenvalues  $\alpha_1=2.38173(1)$ and $\alpha_2=1.60857(1)$ yielding $1/m_{\rm T}=\sqrt{\alpha_1 \alpha_2}=1.95735(1)$.

In a periodic system, the different band edge masses would all be exactly equal $m=m_\rho=m_\text{T}$. Here it is only approximatively true for the clean quasiperiodic case with $\frac{1}{m}\simeq \frac{1}{m_\rho}\simeq \frac{1}{m_\text{T}}\simeq 1.95(2)$ and is not true for the disordered case for which $\frac{1}{m}\simeq \frac{1}{m_\rho}\simeq 1.9 \gg \frac{1}{m_\text{T}} \to 0$.

\section*{Anisotropy of the effective mass tensor}
As we now show using a variational approach, the anisotropy in the mass tensor \mbox{($\alpha_1/ \alpha_2 \simeq 1.481$)}  can be understood as resulting from the anisotropy in link orientations in the tilings. 
{\bf An effective low-energy and long-wavelength dispersion can be obtained from the average slope in the cut and project method. We find that this dispersion relation 
\be
E(\boldsymbol{k})\simeq -2\sum_{j=1}^3 t_j \cos(\boldsymbol{k} \cdot \boldsymbol{\delta}_j),
\ee
is that of an anisotropic triangular lattice with hopping amplitudes $t_j$ (where the nearest neighbor vectors are defined as $\boldsymbol{\delta}_1=(1,0)$, $\boldsymbol{\delta}_2=(-\frac{1}{2},\frac{\sqrt{3}}{2})$ and $\boldsymbol{\delta}_3=(-\frac{1}{2},-\frac{\sqrt{3}}{2})$) reflecting the stoichiometry of the link orientations: $t_1=1-\theta^{-1}$, $t_2=1-\theta^{-3}$, $t_3=1-\theta^{-2}$. In the long wavelength limit $|\boldsymbol{k}|\ll 1/\sqrt{N}$, the dispersion relation can be approximated by a parabola $E(\boldsymbol{k})\simeq E(0)+\frac{1}{2}\alpha_{ij} k_i k_j$, where $E(0)=-\sum_j t_j=-4$, which defines an inverse effective mass tensor $\alpha_{ij}$. Its two eigenvalues $\alpha_1$, $\alpha_2$ give an average inverse effective mass $\frac{1}{m_\text{cl}}=\sqrt{\alpha_1 \alpha_2}\simeq 1.9715$ and an anisotropy $\frac{\alpha_2}{\alpha_1}\simeq1.405$. Both quantities are in fair agreement with those derived from the first mini-band (see the main text). It is remarkable that the same long wavelength approximation describes both the clean and the disordered Rauzy tiling as they share the same average slope in the cut and project method.

The LL width $\Delta E_n$ also contains interesting information. The relative width $\Delta E_n/(E_{n+1}-E_n)$ goes to zero when $f\to 0$, so that LL are almost as well defined as in a periodic crystal for which it is known that the width $\propto \exp(-\# /f)$ (see \cite{Claro79} and references therein). Also, there are obvious cusps in the relative width that correspond to particular fluxes such as $f=\theta^{-2}$. 

Also the LL width has increased. In particular, the relative width now diverges in the zero flux limit, as expected for continuum LL that are broadened by disorder ($\Delta E_n \propto \sqrt{f}$ \cite{Ando74}). These are the Landau bands familiar in the context of IQHE in a weakly disordered 2D electron gas. Cusps in the relative width have also been washed out by phason disorder.

}
\end{document} 
