Efficient algorithm to compute the Berry conductivity

Topological insulators (TI) are a topological state of quantum matter that constitutes a new paradigm in condensed matter physics [1–4]. These recently discovered new materials exhibit unique, fascinating properties such as current-carrying surfaces and edge states, which are strongly protected against perturbations in either the bulk or the surface of the material [5–10], and non-standard exchange statistics of quasi-particle excitations, which offer potential applications in the context of quantum computation [11–13].

The question of what happens in topological insulators when the Fermi energy no longer lies inside the gap between two energy bands is by no means rhetorical but of high practical importance. In fact, this situation naturally occurs in the experimental process of production of candidate samples of topological insulators such as Bi₂Se₃ and Bi₂Te₃ compounds. These are used, for instance, in cooling devices due to their favorable thermoelectric properties. The chemical composition can be well-controlled and adjusted to the composition of the desired topological insulator. However, it is much more demanding to control the level of the Fermi energy, which for many samples lies within the bulk energy bands instead of the insulating energy gap, thereby invalidating them as true TIs. This difficulty has motivated the development of sophisticated molecular beam epitaxy (MBE) techniques to precisely control the growth of ultra-thin Bi₂Se₃ and Bi₂Te₃ films [14, 15]. Likewise, in two-dimensional TIs it is possible to adjust the Fermi energy to lie either in the band gap or the bulk bands. Experimentally, in CdTe/HgTe/CdTe quantum wells, formed by a thin layer of HgTe embedded between two CdTe layers, this can be achieved by an elaborated MBE technique that allows one to control the thickness of the intermediate HgTe layer and thereby tune the position of the Fermi energy with respect to the bands [16, 17]. For an appropriate thickness, the Fermi energy lies in the gap between the bulk bands, and the heterostructure shows the desired characteristic topological insulating behavior with a quantized spin conductivity of $2{{e}^{2}}/h$.

Complementary to solid-state realizations, cold atoms in optical lattices have been proposed as a realistic platform to experimentally explore the new physics of TIs under controllable conditions [18–34]. In particular, in these systems the Fermi energy can be controlled directly by the filling of atoms in the lattice. There are several proposals to measure the transverse conductivity for both the insulating and the metallic case [20, 31, 34]. In contrast, in condensed matter systems such as the above-mentioned chemical compounds, the pinning of the Fermi level to a value inside the bulk bands typically arises due to external causes like crystal defects and other sources, which are not straightforward to control. As a consequence, in transport properties and measurements, bulk carriers often dominate over the contribution stemming from surface or edge states.

Finally, this question also plays a fundamental role in the physics of the anomalous quantum Hall effect (AHE) [35, 36], which precedes the upsurge of topological insulators as a prominent field in condensed matter. In the standard quantum Hall effect (QHE), which can be observed in non-magnetic materials, there is a linear dependence of the Hall resitivity ρ_xy on an externally applied perpendicular magnetic field. In contrast, in the AHE, an anomalous deviation from the linear law is observed in ferromagnetic materials. A complete theory for the AHE has remained elusive for more than a century, largely due to complications arising from the fact that there are three main mechanisms that influence the electronic motion and can give rise to an AHE: the intrinsic mechanism, the skew-scattering mechanism and the side-jump mechanism [36]. Here, we shall be interested in the so-called intrinsic mechanism for the AHE, which is the contribution that can be expressed in terms of the Berry-phase curvature and thereby represents an intrinsic quantum mechanical property of a perfect crystal. This intrinsic contribution, which is dominant in metallic ferromagnets with moderate conductivity, depends only on band structure properties and is largely independent of the scattering that affects other AHE mechanisms.

Understanding this intrinsic and anomalous contribution has become possible with the seminal work by Haldane [37], who uncovered by a fully quantum-mechanical treatment—unlike precedent work based on semiclassical methods [38]—the topological origin of this contribution and its relation to the physics taking place at the Fermi surface. Haldane showed that the intrinsic contribution to the AHE conductivity stems from a combination of an integer-valued part, stemming from the contribution of filled bands, and a part originating from the Fermi surface; i.e., from the cuts of a partially filled band at the Fermi energy E_F (non-integer valued contribution).

It is crucial to realize that in order to directly apply Haldaneʼs equations [37] to a given problem, one needs to know precisely the form of the Fermi surface. In practice, except in very simple model cases, this is not possible, since the band structure of real materials is obtained from detailed numerical calculations. One is typically given a numerical data set about the bands instead of an explicit formula. Thus, in practice it is highly desirable to have at oneʼs disposal a numerical method that is: (i) gauge-invariant, (ii) efficient and (iii) outputs numerical results with controllable and known error intervals. In this work, we develop such a general and efficient numerical algorithm to compute the Berry conductivity when the Fermi energy does not lie within the band gap. In the following, we shall refer to Berry conductivity as the non-quantized conductivity associated with the Chern number according to the Thouless, Kohmoto, Nightingale, Den Nijs (TKNN) formula [39] when the position of the Fermi level lies in the conduction band, so that we recover the TKNN quantized conductivity for the standard insulating case if the Fermi energy lies in the energy gap between two bands.

Our main results are:

• (i)  
  We present a new method to compute the Berry conductivity when the Fermi energy level is located outside the band gap. We outline the algorithm (schematically summarized in figure 1), discuss its ingredients and show that it is gauge invariant and efficient (section 2).
• (ii)  
  We emphasize that a central feature of the presented method is that it is endowed with known and controllable error bars for the non-integer value of the conductivity. This is essential. When the Berry conductivity is not integer-valued, errors due to approximations need to be under control in order to distinguish two different values of the observable conductivity so that one may safely distinguish a topological phase from a trivial phase.
• (iii)  
  To test and benchmark the performance of the algorithm, we first apply it to the paradigmatic Haldane model [40], which has a simple enough structure so that the analytic form of the two-band energy spectrum is known (section 3.1). Subsequently, we apply the method to the more complex case of the Hofstadter model [41], which belongs to the class of multi-band topological insulators, where the band structure information is obtained numerically (section 3.2). These models are both of importance and have attracted interest in the field of quantum simulation of topological insulators with cold atoms in optical lattices. Here, our method provides the theoretical tools that allow one to map out the phase diagrams in future experiments. Finally, we also apply our method to the Bernevig, Hughes and Zhang (BHZ) model [16], which is a realistic model that captures the physics of the 2D spin Hall effect present in systems such as the above-mentioned CdTe/HgTe/CdTe quantum well compounds (section 3.3). We conclude with a short summary and a discussion of possible future extensions of the presented method (section 4).

Zoom In Zoom Out Reset image size

2.1. Generalized Berry conductivity

Before presenting our numerical algorithm to calculate the Berry conductivity, in this section we briefly review the expressions for the intrinsic Hall conductivity for the insulating case where the value of the Fermi energy lies in the gap between two bands. We also review the generalized result for the situation in which the Fermi energy lies in a partially filled band [37].

In the insulating case, the Hall conductivity is quantized and proportional to the sum of the Chern numbers of the occupied energy bands,

$$\begin{eqnarray}&&{{\sigma }_{H}}=\frac{{{e}^{2}}}{h}\mathop{\sum }\limits_{{{E}_{\alpha }}<{{E}_{F}}}{{C}_{\alpha }}.\end{eqnarray} \tag{ 1 }$$

The Chern numbers ${{C}_{\alpha }}$ are integer-valued topological invariants, defined in terms of the integral of the Berry curvature $F_{xy}^{\alpha }({\bf k})$ over the whole Brillouin zone (B.Z.) [39, 42]:

$$\begin{eqnarray}&&\begin{array}{l} {{C}_{\alpha }}=\frac{1}{2\pi i}{{\int }_{B.Z.}}F_{xy}^{\alpha }({\bf k}){{d}^{2}}k \\ F_{xy}^{\alpha }({\bf k})={{\partial }_{{{k}_{x}}}}A_{y}^{\alpha }({\bf k})-{{\partial }_{{{k}_{y}}}}A_{x}^{\alpha }({\bf k}). \\ \end{array}\end{eqnarray} \tag{ 2 }$$

The latter is expressed by the exterior derivative of the Berry connection:

$$\begin{eqnarray}&&A_{\mu }^{\alpha }({\bf k})=\langle {{u}_{\alpha }}({\bf k})|{{\partial }_{\mu }}|{{u}_{\alpha }}({\bf k})\rangle ,\end{eqnarray} \tag{ 3 }$$

where ${{u}_{\alpha }}({\bf k})$ is the eigenvector corresponding to the energy band ${{E}_{\alpha }}({\bf k})$.

In the case in which the Fermi energy does not lie in an energy gap between bands, as schematically shown in figure 1(a), the intrinsic Hall conductivity generalizes to [36, 43]:

$$\begin{eqnarray}&&{{\sigma }_{H}}({{E}_{F}})=\frac{{{e}^{2}}}{h}\mathop{\sum }\limits_{\alpha }{{\mathfrak{C}}_{\alpha }},\end{eqnarray} \tag{ 4 }$$

with

$$\begin{eqnarray}&&{{\mathfrak{C}}_{\alpha }}({{E}_{F}})=\frac{1}{2\pi i}{{\int }_{B.Z.}}{{d}^{2}}k\;F_{xy}^{\alpha }({\bf k})\Theta ({{E}_{F}}-{{E}_{\alpha }}({\bf k})),\end{eqnarray} \tag{ 5 }$$

where $\Theta (E)$ denotes the Heaviside function and α denotes the band index. Thus, the conductivity is the sum of the integer-valued Chern numbers corresponding to fully-occupied energy bands below the Fermi energy as well as a non-quantized contribution that depends on the Fermi surface; i.e., it stems from the integral over energy band(s), which are partially filled at a given Fermi energy E_F .

For systems with a particularly simple band structure (e.g., two-band systems), the expressions for the eigenvalues and eigenvectors of the bands are given in explicit form; hence the Chern values ${{\mathfrak{C}}_{\alpha }}$ can be calculated analytically. In general, however, the Hamiltonian system cannot be diagonalized analytically, and an efficient numerical method to compute the Chern values is needed.

2.2. Construction and properties of the algorithm

The algorithm we propose to numerically compute the Chern values of equation (12), and thereby the Berry conductivity of equation (4), is based on a series of controlled approximations: first, we discretize the two-dimensional Brillouin zone by a finite ${{n}_{B}}\times {{n}_{B}}$ grid of small plaquettes at discrete momenta ${{{\bf k}}_{l}}$ (see figure 1(b) and appendix A for details), so that the integral over the (partially filled) band becomes:

$$\begin{eqnarray}&&{{\mathfrak{C}}_{\alpha }}({{E}_{F}})\longrightarrow \frac{1}{2\pi i}\mathop{\sum }\limits_{\{{{{\bf k}}_{l}}\}}F_{xy,l}^{\alpha }\;p_{l}^{\alpha }({{E}_{F}}),\end{eqnarray} \tag{ 6 }$$

with the Berry curvature contribution:

$$\begin{eqnarray}&&F_{xy,l}^{\alpha }={{\int }_{\square }}{{d}^{2}}k\;F_{xy}^{\alpha }({\bf k})\end{eqnarray} \tag{ 7 }$$

from a small two-dimensional plaquette of size ${{\Delta }_{{{k}_{x}}}}{{\Delta }_{{{k}_{y}}}}$, and the weighting factors:

$$\begin{eqnarray}&&p_{l}^{\alpha }({{E}_{F}})=\frac{1}{{{\Delta }_{{{k}_{x}}}}{{\Delta }_{{{k}_{y}}}}}{{\int }_{\square }}{{d}^{2}}k\;\Theta ({{E}_{F}}-{{E}_{\alpha }}({\bf k})).\end{eqnarray} \tag{ 8 }$$

The weights $p_{l}^{\alpha }({{E}_{F}})$ correspond to the partial area of the plaquette, which is covered by the Fermi sea; thus $p_{l}^{\alpha }({{E}_{F}})=0$ ($p_{l}^{\alpha }({{E}_{F}})=1$) for squares with energies completely above (below) the Fermi energy E_F , and $0<p_{l}^{\alpha }({{E}_{F}})<1$ for momentum space plaquettes that are cut by the Fermi surface (see figure 1(b)). The choice of the value n_B , i.e., the resolution of the momentum space grid, is important: it can be motivated either by given physical conditions, such as a finite experimental energy resolution (or, e.g., the finite size of real-space optical lattices, which in turn induces a smallest characteristic scale in momentum space); or it can be chosen according to given numerical resources. In appendix F, we derive a convergence criterion in terms of the grid discretization and provide an error bound due to the grid discretization.

The key of the numerical algorithm is now to evaluate reliably and, under controlled approximations, the discretized sum of equation (6), whose value converges to equation (4) for increasingly finer grids.

(i) Gauge-invariant calculation of the Berry curvature. To numerically calculate the Berry curvature contributions $F_{xy,l}^{\alpha }$, we employ a numerical algorithm proposed by Fukui, Hatsugai and Suzuki [44] (FHS algorithm). It is highly efficient, and the discrete sum $1/(2\pi i){{\sum }_{\{{{{\bf k}}_{l}}\}}}\tilde{F}_{xy,l}^{\alpha }$ converges rapidly to the correct integer-valued Chern numbers ${{C}_{\alpha }}$, even for a very coarse-grained discretization of the Brillouin zone. This behavior is rooted in the fact that the algorithm is based on a lattice gauge formulation [45, 46] instead of a finite difference discretization of the Berry curvature. In appendix A we provide a brief summary of the FHS algorithm and the explicit expressions for the lattice strength $\tilde{F}_{xy,l}^{\alpha }$ calculated with the FHS method.

(ii) Efficient estimation of the weights $p_{l}^{\alpha }({{E}_{F}})$. To decide whether a given plaquette in momentum space contributes entirely, partially or not at all, we use a simple and rapid classical Monte-Carlo technique: for each plaquette of the grid localized around the discrete momentum ${{{\bf k}}_{l}}$, we generate n_R uniformly distributed random points ${{{\bf k}}_{R}}$ within the plaquette and compute ${{E}_{\alpha }}({{{\bf k}}_{R}})$, which lies above or below the Fermi energy. Based on the latter, we define the estimator:

$$\begin{eqnarray}&&\hat{p}_{l}^{\alpha }({{E}_{F}})=\frac{1}{{{n}_{R}}}\mathop{\sum }\limits_{\{{{{\bf k}}_{R}}\}}\Theta ({{E}_{F}}-{{E}_{\alpha }}({{{\bf k}}_{R}}))\end{eqnarray} \tag{ 9 }$$

for the weighting factors $p_{l}^{\alpha }({{E}_{F}})$.

(iii) Statistical confidence interval and controlled numerical error of the Berry conductivity. Note that the randomness of this estimation procedure introduces a statistical uncertainty. Note also that the value of the estimators $\hat{p}_{l}^{\alpha }({{E}_{F}})$ is bounded between zero and one. However, it is clear that the statistical error will be largest for partially contributing plaquettes with $\hat{p}_{l}^{\alpha }({{E}_{F}})\sim 1/2$, whereas the uncertainty in $\hat{p}_{l}^{\alpha }({{E}_{F}})$ for plaquettes with energies completely above or completely below the Fermi energy is expected to be much smaller. In order to have a known and minimal statistical error in $\hat{p}_{l}^{\alpha }({{E}_{F}})$, and thus in the Berry conductivity, it is highly desirable that the numerical algorithm takes this effect into account and provides statistical errors that depend on the actual value of the Fermi energy. The quantity $\hat{p}_{l}^{\alpha }({{E}_{F}})$ is the estimator of the fixed, though unknown, parameter p of a binomial distribution $\mathcal{B}({{n}_{R}},p)$, corresponding to the process of tossing n_R times a biased coin. As is discussed in detail in appendix B, using the normal approximation, and for a fixed number of runs n_R and a desired value $\epsilon <1$, this allows one to derive a confidence interval $[p_{l,{\rm min} }^{\alpha },\;p_{l,{\rm max} }^{\alpha }]$ for $\hat{p}_{l}^{\alpha }({{E}_{F}})$, called the Wilson interval [47], with modified boundary conditions. This means that with a probability $1-\epsilon$, the 'true' value $p_{l}^{\alpha }$ lies in this interval. The key point is that the width of this interval depends on the actual value of the estimator $\hat{p}_{l}^{\alpha }$ and is typically significantly smaller than the trivial upper bound of one. After symmetrizing the interval by taking the maximum $\Delta \hat{p}_{l}^{\alpha }({{E}_{F}})={\rm max} (\hat{p}_{l}^{\alpha }-{{p}_{{\rm min} }},{{p}_{{\rm max} }}-\hat{p}_{l}^{\alpha })$, each momentum space plaquette of the grid is associated with a probability value ${{\hat{p}}_{l}}({{E}_{F}})\pm \Delta {{\hat{p}}_{l}}({{E}_{F}})$ with a confidence of at least $1-\epsilon$. Finally, we remark that whereas the discussed statistical method is conceptually simple, intuitive and provides controllable error bars, it could be refined and combined with more sophisticated techniques to evaluate the weighting factors (9) or, equivalently, to determine the equal-energy contours of the bands for a given Fermi energy. In addition, an adaptative version of the statistical algorithm could be put forward in which only the weighting factors of plaquettes with large Berry curvature contributions are evaluated with high statistical accuracy.

As mentioned above, for even moderately fine grids, the FHS algorithm provides essentially exact values for the Berry curvature contributions (see [44] and appendix A). Thus, the statistical uncertainty of $\hat{p}_{l}^{\alpha }$ directly translates into an uncertainty in the Berry conductivity contributions,

$$\begin{eqnarray}&&p_{l}^{\alpha }({{E}_{F}})F_{xy,l}^{\alpha }\;\longrightarrow \;\hat{p}_{l}^{\alpha }({{E}_{F}})\tilde{F}_{xy,l}^{\alpha }\pm (\Delta \hat{p}_{l}^{\alpha }({{E}_{F}}))\tilde{F}_{xy,l}^{\alpha }.\end{eqnarray} \tag{ 10 }$$

Finally, the estimated Berry conductivity is given by:

$$\begin{eqnarray}&&{{\tilde{\sigma }}_{{\rm Be}}}({{E}_{F}})=\frac{{{e}^{2}}}{h}\mathop{\sum }\limits_{\alpha }{{\tilde{\mathfrak{C}}}_{\alpha }}({{E}_{F}}),\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&{{\tilde{\mathfrak{C}}}_{\alpha }}({{E}_{F}})=\frac{1}{2\pi i}\mathop{\sum }\limits_{\{{{{\bf k}}_{l}}\}}\tilde{F}_{xy,l}^{\alpha }\;\hat{p}_{l}^{\alpha }({{E}_{F}})\end{eqnarray} \tag{ 12 }$$

with an error $\pm \Delta {{\tilde{\mathfrak{C}}}_{\alpha }}({{E}_{F}})$ of:

$$\begin{eqnarray}&&\Delta {{\tilde{\mathfrak{C}}}_{\alpha }}({{E}_{F}})=\sqrt{\mathop{\sum }\limits_{\{{{{\bf k}}_{l}}\}}{{(\Delta \hat{p}_{l}^{\alpha }({{E}_{F}})\tilde{F}_{xy,l}^{\alpha })}^{2}}}\end{eqnarray} \tag{ 13 }$$

with confidence of at least $1-\epsilon$. We remark that controllable error bars are particularly important and valuable outside the insulating regime; i.e., where the Fermi energy cuts a partially filled energy band. In this case the Berry conductivity is not quantized and can assume continuous non-integer values.

In this section, we apply the algorithm to different models. We first start with the Haldane model, a two-band model that can realize both topological and trivial phases. We then go to the Hofstadter model, a multi-band model characterized by a non-zero Chern number, and finish with the BHZ model, a two-band realistic model realizing a quantum spin Hall effect in condensed matter physics.

3.1. The Haldane model

The model proposed by Haldane [40] is a tight-binding Hamiltonian of spinless fermions on a honeycomb lattice, with dynamics governed by a nearest-neighbor (N.N.) real-valued hopping term of amplitude J and an imaginary next-to-nearest neighbor (N.N.N.) hopping term J₂ (see figure 2(a)). In addition, the fermions are exposed to an onsite staggering potential β, which induces a chemical potential difference between N.N. sites of the bi-partite hexagonal lattice (ϕ and ψ sites). The model is exactly solvable and represents a paradigmatic model in the field of topological phases of matter, as it hosts a quantum AHE phase even in the absence of an external magnetic field. Recently, it has been proposed that the physics of this model could be observed experimentally in a quantum simulation with cold atoms in optical lattices [34].

Zoom In Zoom Out Reset image size

The Hamiltonian of the system is given by:

$$\begin{eqnarray}&&H=-J\mathop{\sum }\limits_{\langle i,j\rangle }c_{i}^{\dagger }{{c}_{j}}+i{{J}_{2}}\mathop{\sum }\limits_{\langle \langle i,j\rangle \rangle }{{\nu }_{ij}}c_{i}^{\dagger }{{c}_{j}}+\beta \mathop{\sum }\limits_{i}{{s}_{i}}c_{i}^{\dagger }{{c}_{i}}.\end{eqnarray} \tag{ 14 }$$

Here, $c_{i}^{\dagger }$ and c_i are fermionic creation and destruction operators, ${{\nu }_{ij}}={\rm sgn} [{{({{{\bf d}}_{1}}\times {{{\bf d}}_{2}})}_{z}}]$ and ${{s}_{\phi ,\psi }}=\pm 1$. The vectors ${{{\bf d}}_{1}}$ and ${{{\bf d}}_{2}}$ are oriented along the bonds of the hexagonal unit cell, as shown in figure 2(a). The model can be readily solved by rewriting the Hamiltonian in terms of two-site basis cells $(\phi ,\psi )$ (see, e.g., [48]) such that the hexagonal lattice becomes a triangular lattice of $(\phi ,\psi )$ cells. In the Fourier space, the Hamiltonian is then given by [40]:

$$\begin{eqnarray}H=\mathop{\sum }\limits_{{\bf k}\in B.Z.}{{\hat{\Psi }}^{\dagger }}({\bf k})\left( \begin{array}{ccccccccccccccc} \beta -2{{J}_{2}}f({\bf k}) & -{{A}^{*}}({\bf k}) \\ A({\bf k}) & -\beta +2{{J}_{2}}f({\bf k}) \\ \end{array} \right)\hat{\Psi }({\bf k}).\end{eqnarray} \tag{ 15 }$$

Here, ${{\hat{\Psi }}^{\dagger }}({\bf k})=(c_{\phi }^{\dagger }({\bf k}),c_{\psi }^{\dagger }({\bf k}))$, $A({\bf k})={\rm exp} (i{\bf k}\cdot {{\delta }_{1}})+{\rm exp} (i{\bf k}\cdot {{\delta }_{2}})+{\rm exp} (i{\bf k}\cdot {{\delta }_{3}})$ is expressed in terms of the vectors between nearest neighbors δ₁, δ₂ and δ₃, and $f({\bf k})={\rm sin} [{{{\bf a}}_{1}}\cdot {\bf k}]+{\rm sin} [{{{\bf a}}_{3}}\cdot {\bf k}]+{\rm sin} [({{{\bf a}}_{1}}+{{{\bf a}}_{2}})\cdot {\bf k}]$ is expressed in terms of the lattice vectors ${{{\bf a}}_{1}}$ and ${{{\bf a}}_{2}}$ as shown in figure 2 and defined in appendix C.

Diagonalization of the Hamiltonian readily yields the two-band energy spectrum:

$$\begin{eqnarray}&&{{E}_{\pm }}({\bf k})=\pm \sqrt{|A({\bf k}){{|}^{2}}+{{(\beta -2{{J}_{2}}f({\bf k}))}^{2}}},\end{eqnarray} \tag{ 16 }$$

which is shown in figure 2(b).

For $\beta ={{J}_{2}}=0$, the Hamiltonian corresponds to pure N.N. hopping of fermions, with the characteristic spectrum exhibiting the two inequivalent Dirac cones [49, 50]. A non-zero staggering potential $\beta \ne 0$ induces an imbalance of the fermion density on the ϕ and ψ lattice sites. The formation of a charge-density-wave phase is associated with the opening of a topologically trivial insulating gap in the spectrum. On the other hand, a strong enough N.N.N. hopping term J₂ opens a topologically non-trivial energy gap that signals the transition of the system into an AHE phase characterized by a non-zero Chern number. The size of the energy gap is determined by the formula $\Delta =2|\beta -3\sqrt{3}{{J}_{2}}|$, and for $|\beta |<3\sqrt{3}|{{J}_{2}}|$ the system is in the topological phase.

We will now illustrate the working principle of our algorithm by applying it step by step—as schematically summarized in figure 1(c)—to the Haldane model. To this end, we start by fixing the Hamiltonian parameters to ${{J}_{2}}=0.1J$, $\beta =0$; i.e., deep in the topologically non-trivial phase. Next, we discretize the Brillouin zone (step 1), where we use for numerical convenience a rectangular-shaped B.Z. parametrization, which is equivalent to the standard hexagonal form (see appendix C for details).

Then, we compute the field strength $\tilde{F}_{xy}^{\alpha }$ for each plaquette (step 2); the result is shown in the right column of figure 3. We fix the number of random points (we choose n_R = 20) (step 3) and compute for each plaquette for n_R randomly distributed momentum vectors ${{E}_{\alpha }}({{{\bf k}}_{R}})$ (step 4). Once the Fermi energy is fixed (step 5), here to a value of ${{E}_{F}}=-1.5J$ so that the Fermi energy level cuts the lower band, we compute the estimators for the weights $\hat{p}_{l}^{\alpha }({{E}_{F}})$ according to equation (9) (step 6). The values of the estimators $\hat{p}_{l}^{\alpha }({{E}_{F}})$ are shown in the left column of figure 3. The central column of the figure displays the associated statistical uncertainties $\Delta \hat{p}_{l}^{\alpha }({{E}_{F}})$, as determined in step 7 with the Wilson interval with modified boundaries, and symmetrized (see appendix B). As expected and desired, the statistical errors associated with plaquettes that correspond to regions that clearly lie above or below the Fermi energy are minimal. In constrast, the plaquettes at energies around E_F , which are cut by the Fermi surface, have higher values. Note that even for very limited Monte Carlo statistics involving only n_B = 20 random points per momentum space plaquette, these uncertainty values are still much smaller than the upper bound of one. In fact, higher uncertainties and error bars for plaquettes around the Fermi surface reflect the physical fact that these are the plaquettes corresponding to the regions in momentum space where small changes in the Fermi energy level can lead to smaller or larger contributions of Berry curvature and thus to changes in the Berry conductivity. The central and lower row of figure 3 show the weight estimators, uncertainties and Berry curvature contributions for larger values of n_B , illustrating how an increasingly finer grid of the Brillouin zone leads to increased resolution and numerical precision.

Zoom In Zoom Out Reset image size

Finally, the estimated weights $\hat{p}_{l}^{\alpha }({{E}_{F}})$ and the Berry curvature contributions $\tilde{F}_{xy,l}^{\alpha }$ are combined to calculate the Berry conductivity (step 8) according to equations (11) and (12), with an associated error bar (step 9) as given by equation (13). By applying the algorithm again for varying values of the Fermi energy, the Berry conductivity can be obtained as a function of the Fermi level energy. The obtained Berry conductivity is shown in figure 4(a): starting from low conductivity values at the bottom of the lower energy band, the conductivity increases up to its plateau value of one (for Fermi energies lying in the topological insulating gap) before it subsequently starts to fall off again once the Fermi energy reaches the upper band.

Zoom In Zoom Out Reset image size

To test the behavior of the algorithm when the system undergoes a phase transition from the topological AHE phase to the trivial insulating phase, we increase the Hamiltonian parameter β to observe the competition of the N.N.N. hopping term with the staggering potential. The subplots in figure 4 show the transition from the topologically non-trivial phase characterized by a Chern number of one to the topologically trivial charge-density-wave phase with a vanishing Chern number. The algorithm correctly captures the closing of the gap as well as the jump of the conductivity plateau-value as the phase transition takes place. We emphasize that the algorithm automatically takes into account the fact that, at the phase transition, the Berry curvature is highly localized at the Dirac points and thus concentrated in only a few plaquettes—a fact that the algorithm signals in the form of larger error bars of the Berry conductivity in the parameter regime where the transition occurs.

Finally, we apply the algorithm to the case where the system resides in the topological phase with a small topological gap opened. Here, the algorithm allows one to clearly verify numerically the $1/{{E}_{F}}$ power law dependence of the Berry conductivity for Fermi energies close to the gap. The ${{\sigma }_{Be}}({{E}_{F}})=({{e}^{2}}/h)\;3\sqrt{3}{{J}_{2}}/|{{E}_{F}}|$ behavior is predicted by the linear approximation of the spectrum around the Dirac points [43, 51–53]. The results are shown and discussed in figure 5.

Zoom In Zoom Out Reset image size

3.2. The Hofstadter model

Let us now apply the numerical algorithm to the Hofstadter model [41], which describes spinless fermions on a square lattice subjected to a uniform magnetic field of magnetic flux quanta per unit cell $\Phi$. Only very recently, several groups have achieved the observation of the characteristic physics, including the fractal spectrum known as Hofstadterʼs butterfly in graphene superlattice systems [54–56]. This work complements ongoing experimental efforts to realize theoretical ideas [57] on how to implement the fermionic Hofstadter Hamiltonian with cold atoms in optical lattices [58–61].

The Hamiltonian in second-quantized form is given by:

$$\begin{eqnarray}&&H=\mathop{\sum }\limits_{\langle i,j\rangle }{{e}^{i{{\theta }_{ij}}}}c_{i}^{\dagger }{{c}_{j}},\end{eqnarray} \tag{ 17 }$$

where the sum is over N.N. sites (see figure 6), and the phase factor ${\rm exp} (i{{\theta }_{ij}})$ corresponds to the Peierls substitution, expressed in terms of the line integral over the vector potential along the link between two neighboring sites i and j of the square lattice. If $\Phi =p/q$ is a rational number, the energy spectrum of the bulk, described in the Fourier space, splits into q sub-bands, each one of them associated with a non-trivial integer-valued Chern number.

Zoom In Zoom Out Reset image size

Due to its multi-band structure, the Hofstadter Hamiltonian can in general not be diagonalized analytically and thus represents an interesting testbed for the numerical algorithm. Figure 7 shows the numerical results for the Berry conductivity for different values of the flux per plaquette ($\Phi {\mkern 1mu} =$ 1/3, 1/5 and 1/7). For Fermi energies lying in the energy gap between bulk bands, the algorithm correctly reproduces the constant Berry conductivity, which corresponds to the sum of the Chern numbers of completely filled bands. Once the Fermi energy falls into a bulk band, the Berry conductivity is no longer quantized. Whereas for $\Phi =1/3$ the Berry conductivity interpolates monotonically between the gap plateau values, for $\Phi =1/5$ the conductivity displays an interesting feature for Fermi energy values in the second band: instead of showing monotonic growth, it first decreases to a minimum value before starting to increase until it reaches the plateau dictated by the quantized value of the conductivity in the gap. The same phenomenon occurs, even more markedly, in the third band of the spectrum for $\Phi =1/7$. The small controlled statistical error bars of the numerical method ensure that the non-monotonic signature in the Berry conductivity is indeed a physical feature rather than a numerical artifact.

Zoom In Zoom Out Reset image size

3.3. The BHZ model

In 2005, it was suggested that the quantum spin Hall effect (QSHE) could possibly be observed in graphene [51, 62]. The QSHE, however, turned out to be impeded by too weak spin–orbit coupling in this system. Shortly later, a realization of the QSHE in HgTe/CdTe nanowell structures was proposed [16] and experimentally realized only one year later [17]: by varying the thickness of the different layers of the heterostructure, the material can exhibit a trivial insulating phase as well as a topological insulating phase, characterized by a ${{\mathbb{Z}}_{2}}$ topological invariant. The physics can be described by an effective Hamiltonian valid close to the $\Gamma$ point, derived by Bernevig, Hughes and Zhang (BHZ model) [16, 63]. The Hamiltonian is given by a 4 × 4 matrix in momentum space,

$$\begin{eqnarray}\begin{array}{l} H=\left( \begin{array}{ccccccccccccccc} h({\bf k}) & 0 \\ 0 & {{h}^{*}}(-{\bf k}) \\ \end{array} \right), \\ h({\bf k})=\epsilon ({\bf k})+{{d}_{i}}{{\sigma }^{i}}, \\ \end{array}\end{eqnarray} \tag{ 18 }$$

where $\mathcal{1}$ is the two-dimensional identity matrix, σⁱ denote the Pauli matrices and:

$$\begin{eqnarray}&&\begin{array}{l} \epsilon ({\bf k})=C-D(k_{x}^{2}+k_{y}^{2}), \\ {\bf d}({\bf k})=(A{{k}_{x}},-A{{k}_{y}},M({\bf k})), \\ M({\bf k})=M-B(k_{x}^{2}+k_{y}^{2}). \\ \end{array}\end{eqnarray} \tag{ 19 }$$

The parameters A, B, C, D and M depend on material properties as well as the thickness of the layers and can be computed numerically [16, 63].

The Hamiltonian decouples into 2 × 2 blocks, and the spin conductivity can be written as the difference of the conductivity for each spin orientation, thus it makes sense to study the conductivity of one of the orientations. Here, we apply our algorithm to the BHZ model with parameters as calculated in [16]. Figure 8(a) shows the energy spectrum that exhibits a small gap of $0.01\;{\rm eV}$, which renders the computation of the Berry conductivity in the non-insulating regime more demanding. Figures 8(b)–(d) show the numerical results for the Berry conductivity for increasingly finer grids of the Brillouin zone.

Zoom In Zoom Out Reset image size

Even for the roughest grid studied (${{n}_{B}}=40$), the algorithm correctly captures the qualitative behavior and the conductivity minimum value of $-1$ for the Fermi energy lying in the shallow energy gap. However, as signaled by considerably large error bars, only a few plaquettes contribute large values of Berry curvature to the conductivity. Thus, finer grids are required to quantitatively correctly describe the conductivity behavior in the vicinity of the gap (see figures 8(c) and (d) with n_B = 160 and n_B = 320). This effect illustrates the importance of a high enough resolution, both numerically and in an experiment. As the algorithm qualitatively captures the behavior even for rather coarse-grained grids, this can be helpful to predict observations in the case of restricted experimental resolution, e.g., originating from the finite size of optical lattices for cold atoms or finite temperature constraints in solid state experiments. In appendix F, error bounds for the conductivity that take into account a finite grid resolution are discussed in detail.

Finally, we remark that the BHZ model is an effective model valid close to the $\Gamma$ point, and thus the results of our analysis are also only valid in the vicinity of the energy gap. It is possible and will be an interesting extension of the present work to apply the numerical method to a more realistic, refined model, which incorporates more information about the band structure of the system.

In this work we have proposed and constructed a numerical algorithm to calculate the Berry conductivity in topological band insulators. The algorithm works for the insulating case where the Fermi energy lies in the gap between two bulk bands; it also works for the situation in which it lies within a band. The algorithm is gauge-invariant by construction, efficient and outputs the Berry conductivity with known and controllable error bars. We have successfully applied the algorithm to several paradigmatic models of topological quantum matter, including Haldaneʼs model on the honeycomb lattice [40], the multi-band Hofstadter model [41] and the BHZ model [16] that describes the 2D spin Hall effect observed in CdTe/HgTe/CdTe quantum well compounds.

In addition to its applicability to topological insulators, the numerical method to compute the Berry conductivity for arbitrary values of the Fermi energy level can be applied to several other important problems: it can be used to study new phases of matter such as topological Fermi liquids [37, 64, 65] that arise in interacting systems of fermions that realize a TI phase or an AHE phase. Mean field methods applied to these systems predict the existence of such phases [48, 66, 67]. Here, the efficient and controllable numerical method for computing the Berry conductivity provides the appropriate observable to map out the possible topological phases of those systems with the desired accuracy [68–70].

Recent experiments, in which insulating phases [58, 61] have been quantum simulated with cold atoms in optical lattices, provide another natural scenario where our new algorithm can be applied. Complementary to condensed matter systems, these experimental setups offer the possibility to study the intrinsic Berry conductivity in AHE systems under particularly clean and controllable conditions. Here, our algorithm can provide a precise observable to reliably and quantitatively distinguish symmetry protected topological phases from trivial phases. It also can predict some interesting features within the energy band. In fact, several ways have been proposed to measure characteristic signatures of topological quantum phases in systems of cold atoms [71–75]. In particular, recently several ways to measure the Berry conductivity in cold atom experiments using time of flight measurements have been proposed [20, 31, 34].

An experimentally useful extension of our work would be to generalize our numerical method to the case of three dimensional topological insulators under time-reversal symmetry protecting conditions. Finally, an interesting question is how to generalize the controlled numerical method to an open quantum system scenario, such that it can be applied to topological insulators and topologically ordered systems coupled to an environment [76–80].

AD thanks the FRS-FNRS Belgium for financial support and N Goldman and P Gaspard for support and valuable discussions. We acknowledge support by the Spanish MICINN grants FIS2009-10061 and FIS2012-33152, the CAM research consortium QUITEMAD S2009-ESP-1594, the European Commission PICC: FP7 2007-2013, Grant No. 249958, and the UCM-BS grant GICC-910758.

The continuous Brillouin zone is discretized by a two-dimensional lattice grid of n_B points in each direction. For simplicity, we focus here on a rectangular grid, but the formalism can be readily extended to any polygonal grid [46]. The plaquettes of the momentum space lattice are then given by:

$$\begin{eqnarray}&&{{{\bf k}}_{l}}={{{\bf k}}_{{\rm min} }}+i{{{\bf s}}_{{{k}_{x}}}}+j{{{\bf s}}_{{{k}_{y}}}},\end{eqnarray} \tag{ A.1 }$$

with

$$\begin{eqnarray}&&\begin{array}{l} 0\leqslant i,j\leqslant {{n}_{B}}-1, \\ \quad {{{\bf s}}_{{{k}_{x}}}}=\delta {{k}_{x}}{{{\bf u}}_{{{k}_{x}}}}, \\ \quad {{{\bf s}}_{{{k}_{y}}}}={{\delta }_{{{k}_{y}}}}{{{\bf u}}_{{{k}_{y}}}}, \\ \quad {{\delta }_{{{k}_{x}}}}=({{k}_{xmax}}-{{k}_{xmin}})/{{n}_{B}}, \\ \quad {{\delta }_{{{k}_{y}}}}=({{k}_{ymax}}-{{k}_{ymin}})/{{n}_{B}}. \\ \end{array}\end{eqnarray} \tag{ A.2 }$$

The lattice field strength $\tilde{F}_{xy}^{\alpha }({{{\bf k}}_{l}})$ of band α on the grid is then defined in terms of the link variable ${{U}_{\mu }}({\bf k})$ as:

$$\begin{eqnarray}&&{{\tilde{F}}_{xy}}({{{\bf k}}_{l}}):={\rm ln} [{{U}_{x}}({{{\bf k}}_{l}}){{U}_{y}}({{{\bf k}}_{l}}+{{{\bf 1}}_{{{k}_{x}}}})/({{U}_{x}}({{{\bf k}}_{l}}+{{{\bf 1}}_{{{k}_{y}}}}){{U}_{{{k}_{y}}}}({{{\bf k}}_{l}}))],\end{eqnarray} \tag{ A.3 }$$

where ${{U}_{\mu }}({\bf k})=\langle u({\bf k})|u({\bf k}+{{{\bf 1}}_{\mu }})\rangle$. If the admissibility condition $|{{\tilde{F}}_{xy}}({{{\bf k}}_{l}})|<\pi$ is satisfied [44, 46], the lattice gauge theory corresponds to the continuous gauge theory [44, 46], and one can write:

$$\begin{eqnarray}&&|{{F}_{xy}}({{{\bf k}}_{l}})|\delta {{k}_{x}}\delta {{k}_{y}}\simeq |{{\tilde{F}}_{xy}}({{{\bf k}}_{l}})|.\end{eqnarray} \tag{ A.4 }$$

Based on these Berry curvature contributions, the Chern number can be computed as:

$$\begin{eqnarray}&&\tilde{C}=\frac{1}{2\pi i}\mathop{\sum }\limits_{{{{\bf k}}_{l}}}{{\tilde{F}}_{xy}}({{{\bf k}}_{l}}).\end{eqnarray} \tag{ A.5 }$$

In this section, we present the concept and the details of a confidence interval (C.I.) to characterize the statistical uncertainty of the estimated weights $\hat{p}_{l}^{\alpha }$, as defined in equation (9). For simplicity of the notation, we suppress the band index α and momentum index l in the following.

The problem of estimating the weights corresponds to determining the unknown, though fixed probability value p of a binomial distribution $\mathcal{B}({{n}_{R}},\;p)$, based on the outcome of n_R trials. The probability to observe k of the n_R enquiries the value +1 is given by:

$$\begin{eqnarray}&&P(X=k)=\frac{{{n}_{R}}!}{k!({{n}_{R}}-k)!}{{p}^{k}}{{(1-p)}^{k}}.\end{eqnarray} \tag{ B.1 }$$

The goal is to associate a C.I. of a width much smaller than one with the estimated value $\hat{p}$, such that the true value p lies with a probability $1-\epsilon$ inside the C.I. There are several ways to define the C.I., and we will in the following outline the advantages and inconveniences of some of them to motivate the necessity to adopt a simple and appropriate one that we use in our algorithm. To characterize and compare the quality of different conventions for the C.I, it is convenient to introduce the coverage probability: it corresponds to the effective probability to be inside the C.I. and can be compared to the expected probability $1-\epsilon$. As a guiding principle, a 'good' C.I. is an interval with ${{p}_{\operatorname{cov}}}\simeq 1-\epsilon$. On the contrary, for ${{p}_{\operatorname{cov}}}<1-\epsilon$, the C.I. is 'bad' as the statistical 'guaranteeing functionality' of the interval fails. The other case ${{p}_{\operatorname{cov}}}>1-\epsilon$ is not dramatic in our context, as this implies that the true value of the estimated quantity p actually lies in the C.I. with a probability even higher than the targeted value of $1-\epsilon$.

The construction of the C.I. is based on the central limit theorem, which can be used to prove the convergence of the binomial distribution to a normal distribution $\mathcal{N}$, in our case:

$$\begin{eqnarray}&&\sqrt{{{n}_{R}}}\frac{{{{\hat{p}}}_{l}}-p}{\sqrt{p(1-p)}}\to \mathcal{N}(0,1)\qquad \quad {\rm when}\;{{n}_{R}}\to \infty .\end{eqnarray} \tag{ B.2 }$$

The central limit theorem and the definition of the C.I. of a normal distribution with an expected probability $1-\epsilon$ permits us to write the C.I. of ${{\hat{p}}_{l}}$ as a self-consistent equation in terms of p:

$$\begin{eqnarray}&&p={{\hat{p}}_{l}}\pm {{z}_{\epsilon /2}}\sqrt{\frac{p(1-p)}{{{n}_{R}}}},\end{eqnarray} \tag{ B.3 }$$

where ${{z}_{\epsilon }}$ is the quantile function of the normal distribution [81].

A first way to define a C.I. is by maximizing the second term of the sum, yielding:

$$\begin{eqnarray}&&p={{\hat{p}}_{l}}\pm \frac{{{z}_{\alpha /2}}}{2\sqrt{{{n}_{R}}}}\end{eqnarray} \tag{ B.4 }$$

for $p=1/2$. This relation highlights the typical $1/\sqrt{{{n}_{R}}}$ dependence of the statistical error and can be used to provide a rough estimate of the size of the C.I. in terms of n_R . However, as the length of the interval no longer depends on the estimated value ${{\hat{p}}_{l}}$ itself, it does not satisfy our requirement. It will have a coverage probability ${{p}_{\operatorname{cov}}}>1-\epsilon$ and would output error bars that overestimate the actual uncertainty of the observable of interest.

Another commonly used C.I. is constructed using the approximation $p(1-p)\simeq {{\hat{p}}_{l}}(1-{{\hat{p}}_{l}})$ in equation (B.3), thereby replacing the unknown 'true' value by the estimator value, so that:

$$\begin{eqnarray}&&p={{\hat{p}}_{l}}\pm {{z}_{\epsilon /2}}\sqrt{\frac{{{{\hat{p}}}_{l}}(1-{{{\hat{p}}}_{l}})}{{{n}_{R}}}}.\end{eqnarray} \tag{ B.5 }$$

This C.I. is known as the Wald interval [47]. Despite its simplicity, this convention suffers from several problems: for ${{\hat{p}}_{l}}\simeq 0$ or ${{\hat{p}}_{l}}\simeq 1$, the Wald interval shrinks to zero, implying a bad coverage probability for p-values close to one or zero. As discussed by Brown et al [47], a series of criteria has been used in the literature to test the region of validity of this C.I. However, these criteria can be misleading and do not always characterize correctly the C.I.

In figure B.1(a) we illustrate this problem for a fixed value of n_R = 40 and by computing the coverage probability of the C.I. in terms of the value of p on 10 000 samples. One notices at first glance the tendency of the curve to lie below the expected value of $1-\epsilon$. Although the C.I. works rather well for values of p close to p = 0.5, it captures only poorly the situation at values close to the boundaries. Finally, the curve has fast and significant oscillating behavior, which gives rise to the phenomenon of so-called 'lucky'/'unlucky' numbers when increasing slightly the probability p, the coverage probability jumps from a good ${{p}_{\operatorname{cov}}}$ to a poor ${{p}_{\operatorname{cov}}}$ value, as is the case, for instance, around p = 0.8 in the example shown. The couple $(p,{{n}_{R}})$ defines the lucky/unlucky numbers. In figure B.2(a) we fix the value p = 0.25 and vary the value of n_R . Here one also observes significant fluctuations that are only stabilized at larger values of n_R . This effect becomes even more striking at small p, as illustrated in figure B.2(c), where a fixed value of p = 0.007 has been chosen: under an increase of n_R , the C.I. seems to converge to a favorable value of ${{p}_{\operatorname{cov}}}$ until reaching n_R = 423, where ${{p}_{\operatorname{cov}}}$ suddenly drops from 0.94 to 0.78. We thus exclude the Wald interval as a candidate to construct the C.I. for the $\hat{p}_{l}^{\alpha }$ estimators in our algorithm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Most of the mentioned problems can be avoided if the approximation $p(1-p)\simeq {{\hat{p}}_{l}}(1-{{\hat{p}}_{l}})$ is not applied in equation (B.3). Instead, one can exactly solve equation (B.3), which is a quadratic equation for $\hat{p}$. This yields the so-called Wilson interval [47, 82]:

$$\begin{eqnarray}&&{{p}_{{\rm max} ,{\rm min} }}=\left( {{{\hat{p}}}_{l}}+\frac{z_{\epsilon /2}^{2}}{2{{n}_{R}}}\pm {{z}_{\epsilon /2}}\sqrt{\frac{{{{\hat{p}}}_{l}}(1-{{{\hat{p}}}_{l}})}{{{n}_{R}}}+\frac{z_{\epsilon /2}^{2}}{4n_{R}^{2}}} \right)/\left( 1+\frac{z_{\epsilon /2}^{2}}{{{n}_{R}}} \right).\end{eqnarray} \tag{ B.6 }$$

As illustrated in figure B.1(b), the Wilson interval is much more stable, and the coverage probability is oscillating around the value $1-\epsilon$. Figure B.2(b) shows that the Wilson interval reaches rapidly and in a stable way the expected value $1-\epsilon$. The only problem still to be cured is at the boundaries, at p-values around zero or one, where the coverage probability drops. Figure B.2(d) illustrates the convergence at small p, here fixed to p = 0.007, and indicates that the effect of lucky/unlucky numbers is much less important than for the Wald interval. Brown et al [47] propose to replace the lower (upper) boundary of the C.I., obtained by the normal approximation by a lower (upper) boundary obtained from a Poisson approximation for small (big) values of $\hat{p}$. This indeed stabilizes the behavior of the C.I. even close to the boundaries, but complicates the expression of the C.I. Here, we propose a simpler patch, which has the same desired effect: we use the following replacement:

$$\begin{eqnarray}\begin{array}{rcl} {{p}_{{\rm min} }} & = & 0\quad \;{\rm if}\;{{{\hat{p}}}_{l}}=x/{{n}_{R}},x=0,1,2, \\ {{p}_{{\rm max} }} & = & 1\quad \;{\rm if}\;{{{\hat{p}}}_{l}}=x/{{n}_{R}},x={{n}_{R}},{{n}_{R}}-1,{{n}_{R}}-2, \\ \end{array}\end{eqnarray} \tag{ B.7 }$$

including x = 3 and $x={{n}_{R}}-3$ if ${{n}_{R}}>40$. Finally, merely for convenience to obtain symmetric error bars, we symmetrize the C.I. around $\hat{p}$ by choosing a width that corresponds to twice the value of ${\rm max} \{{{p}_{{\rm max} }}-\hat{p},\hat{p}-{{p}_{{\rm min} }}\}$. While keeping the C.I. narrow, this only leads to a modest over-estimation of the actual uncertainty of the estimator.

The C.I. interval defined in this form has a simple analytical form in combination with a good coverage probability, even for small n_R [47]. We will use this construction of the C.I. in the Monte Carlo sampling part of the algorithm and refer to it as 'Wilson interval with modified boundaries' in the main text.

The reciprocal vectors are ${{{\bf b}}_{1}}=(2\pi /3,2\pi /\sqrt{3})$ and ${{{\bf b}}_{2}}=(-2\pi /3,2\pi /\sqrt{3})$. The equivalence between the hexagonal Brillouin zone and the rectangular area used in the computation is shown in figure C.1: the lower left triangle can be translated along ${{{\bf b}}_{1}}$ and the lower right triangle can be translated along the ${{{\bf b}}_{2}}$.

Zoom In Zoom Out Reset image size

In this section, we illustrate the importance of an appropriate momentum space resolution, parametrized by the discretization number n_B . Figure C.2 presents a zoom of figure 4(c) of the main text, showing the numerically estimated Berry curvature for different grids n_B . One finds that all graphs have the same behavior until reaching a value around ${{E}_{F}}=-0.33J$. There, the behavior of the estimator of the Berry curvature becomes jerky. This is a characteristic which shows up when a few momentum-space plaquettes have an important Berry curvature contribution. The error bars signal this effect. The situation improves for increasing values of n_B : the curves converge to one sharp curve, showing that the main contribution of the Berry curvature stems from states with an energy close to zero, and the error bars decrease significantly.

Zoom In Zoom Out Reset image size

Another way to reduce the size of the error bars is to increase the value of n_R , the number of random points used to compute ${{\hat{p}}_{l}}$ in each plaquette. Figure D.1 displays the estimated Berry conductivity for a fixed value n_B = 20 and for the two values n_R = 20 and n_R = 160. As expected, the curve for n_R = 160 is much more stable than the curve obtained for n_R = 20: we see here a better interpolation in terms of the Fermi energy at this resolution. We emphasize the fact that the curve corresponding to n_R = 160 is contained completely in the region spanned by the error bars of the n_R = 20. This is an important point of the chosen construction of the confidence interval, as described in appendix B.

Zoom In Zoom Out Reset image size

In this section, we compare the Wilson interval with modified boundaries with the C.I. defined in equation (B.4) by examining the final error interval obtained in the Haldane model using both methods. We work here with ${{J}_{2}}=0.1J$, $\beta =0.5J$ such that the Berry curvature is really sharp and localized. Figure E.1 shows the results for both types of error interval with the parameters n_B = 20, n_R = 100. The error interval as obtained by using the Wilson interval (see appendix B) captures correctly the fact that the main contribution to the Berry curvature is strongly localized in momentum space. This gives rise to an increased statistical error in the energy region, in which the Fermi energy crosses plaquettes with a large contribution to the Berry curvature. The error obtained with the other C.I., presented in figure E.1(b), is constant, independent of the value of the Fermi energy, and is thus not indicating the region where the Fermi energy crosses plaquettes with a large contribution to the Berry curvature. This point illustrates the choice of the Wilson interval to construct the main algorithm.

Zoom In Zoom Out Reset image size

The algorithm introduced in this paper allows one to compute the Berry conductivity with controllable statistical error bars for a given grid resolution of the discretization of the Brillouin zone.

The choice of this grid can be dictated either by physical constraints of the problem, such as the finite size of the considered lattices in real space (e.g., in experiments with cold atoms in optical lattices), or by limited numerical resources. In any case, it is important and highly desirable to be able to characterize the error due to the grid resolution. More precisely, for a given grid, the aim is to provide an upper bound on how much the Berry conductivity might at most change if the grid chosen were even finer. In this section, we address this problem and construct a tight error bound in terms of the grid resolution. In particular, this bound will require no inputs except the Berry curvature contributions of the small momentum space plaquettes, which in any case need to be determined (see step 2 in figure 1(c)) in the course of computing the Berry conductivity by our algorithm.

We start with a discretization of the B.Z. into L × L plaquettes such that the sum $\sum _{{{i}_{1}}=1}^{{{L}^{2}}}{{\tilde{F}}_{{{i}_{1}}}}$ over all the plaquettes is equal to the Chern number. (In this section we are omitting the band index for notational simplicity.) Consider then finer and finer grids, where the nth level grid contains ${{c}^{n-1}}L\times {{c}^{n-1}}L$ plaquettes at each iteration, the previous plaquette is split into c² new small plaquettes as illustrated in figure F.1, which presents two successive grid levels with c = 2.

Zoom In Zoom Out Reset image size

The Berry conductivity at the iteration n is written as:

$$\begin{eqnarray}&&{{\tilde{\mathfrak{C}}}^{(n)}}=\sum 1=1{{L}^{2}}i\sum 2=1{{c}^{2}}i\cdots \sum n=1{{c}^{2}}i{{\tilde{F}}_{{{i}_{1}}\ldots {{i}_{n}}}}{{p}_{{{i}_{1}}\ldots {{i}_{n}}}},\end{eqnarray} \tag{ F.1 }$$

where

$$\begin{eqnarray}&&{{\tilde{F}}_{{{i}_{1}}\ldots {{i}_{n}}}}=\frac{1}{{{c}^{2}}}{{\tilde{F}}_{{{i}_{1}}\ldots {{i}_{n-1}}}}(1+{{\epsilon }_{{{i}_{1}}\ldots {{i}_{n}}}})\end{eqnarray} \tag{ F.2 }$$

and

$$\begin{eqnarray}&&\frac{1}{{{c}^{2}}}\sum n=1{{c}^{2}}i{{p}_{{{i}_{1}}\ldots {{i}_{n}}}}={{p}_{{{i}_{1}}\ldots {{i}_{n-1}}}}.\end{eqnarray} \tag{ F.3 }$$

The parameters ${{\epsilon }_{{{i}_{1}}\ldots {{i}_{n}}}}$ quantify the non-homogeneous contribution of the Berry curvature of the subplaquette in terms of the Berry curvature of the plaquette of the previous iteration.

The error bound at the iteration n,

$$\begin{eqnarray}&&\Delta \tilde{\mathfrak{C}}_{\infty }^{(n)}:=\mathop{\sum }\limits_{m=0}^{\infty }\Delta {{\tilde{\mathfrak{C}}}^{(n+m)}},\end{eqnarray} \tag{ F.4 }$$

is defined as the sum over all the relative errors $\Delta {{\tilde{\mathfrak{C}}}^{(n+m)}}={{\tilde{\mathfrak{C}}}^{(n+m)}}-{{\tilde{\mathfrak{C}}}^{(n+m+1)}}$ between the iterations $n+m$ and $n+m+1$.

To compute the upper bound of this quantity, we should make an assumption about the smoothness of the Berry curvature; we expect that when the coarse graining is sufficiently fine after n₀ iterations, the Berry curvature becomes smoother at each iteration. Formally, we assume that there exist an n₀ and a parameter $0<q<1$ such that for $n>{{n}_{0}}$:

$$\begin{eqnarray}&&{{\epsilon }^{(n)}}:=\mathop{{\rm max} }\limits_{{{i}_{1}}\cdots {{i}_{n}}}|{{\epsilon }_{{{i}_{1}}\ldots {{i}_{n}}}}|\leqslant q\mathop{{\rm max} }\limits_{{{i}_{1}}\cdots {{i}_{n-1}}}|{{\epsilon }_{{{i}_{1}}\ldots {{i}_{n-1}}}}|.\end{eqnarray} \tag{ F.5 }$$

Note that it is essential to numerically verify for a given model and set of Hamiltonian parameters that this natural assumption is indeed fulfilled, and to determine from which n₀ on (see also examples below).

Using the last inequality, it is straightforward to derive an upper bound of the Berry conductivity ${{C}^{(n+1)}}$ and of the relative error $\Delta {{C}^{(n)}}$ in terms of ${{C}^{(n)}}$:

$$\begin{eqnarray}&&{{\tilde{\mathfrak{C}}}^{(n+1)}}\leqslant (1+q{{\epsilon }^{(n)}}){{\tilde{\mathfrak{C}}}^{(n)}},\end{eqnarray} \tag{ F.6 }$$

$$\begin{eqnarray}&&\Delta {{\tilde{\mathfrak{C}}}^{(n)}}\leqslant q{{\tilde{\mathfrak{C}}}^{(n)}}.\end{eqnarray} \tag{ F.7 }$$

By iterating, one finds the bound of $\Delta {{\tilde{\mathfrak{C}}}^{(n+m)}}$ in terms of ${{\tilde{\mathfrak{C}}}^{(n)}}$:

$$\begin{eqnarray}&&\Delta {{\tilde{\mathfrak{C}}}^{(n+m)}}\;\leqslant {{q}^{m+1}}{{\epsilon }^{(n)}}{{(1+q{{\epsilon }^{(n)}})}^{m}}{{\tilde{\mathfrak{C}}}^{(n)}}.\end{eqnarray} \tag{ F.8 }$$

The total error $\Delta \tilde{\mathfrak{C}}_{\infty }^{(n)}$ is bounded by a geometric series which, when $q(1+q{{\epsilon }^{(n)}})<1$, converges and can be bounded by the value:

$$\begin{eqnarray}&&\Delta \tilde{\mathfrak{C}}_{\infty }^{(n)}\leqslant \frac{q}{1-q(1+q{{\epsilon }^{(n)}})}{{\epsilon }^{(n)}}{{\tilde{\mathfrak{C}}}^{(n)}}.\end{eqnarray} \tag{ F.9 }$$

As an illustration, we apply this formula to the Haldane model in the case of ${{J}_{2}}/J=0.5,\beta /J=0$. In this case, the Berry curvature is smooth, and the FHS algorithm already converges for n_B = 10. We thus choose $L={{n}_{B}}=10$ for the first iteration. Figure F.2(a) shows the Berry curvature after two iterations with c = 2 (i.e., n_B = 40). Following the exposed line of argument, we first verify numerically that we are in a convergence regime by testing the inequality (F.5) at each iteration. Figures F.2(b) and (c) present the maximum of the ${{\epsilon }_{{{i}_{1}}\ldots {{i}_{n}}}}$ for two successive iterations in terms of the plaquette of the grid n_B = 40. As expected, we find that the contribution of the at the next iteration level is getting smaller. In table 1, we present more quantitative results for the different iterations for the $\epsilon _{{\rm max} }^{(n)}$ at each iteration and the value of ${{q}_{n}}=\epsilon _{{\rm max} }^{(n+1)}/\epsilon _{{\rm max} }^{(n)}$. As can be seen from the table, the value of ${{q}_{n}}=\epsilon _{{\rm max} }^{(n+1)}/\epsilon _{{\rm max} }^{(n)}$ is inferior to 1, and by choosing q = 0.6, we ensure that relation (F.5) is satisfied.

Zoom In Zoom Out Reset image size

Table 1.  The convergence of the Haldane model is studied for the parameter $\beta /J=0$ and two values of the parameter J₂: ${{J}_{2}}/J=0.5$ and ${{J}_{2}}/J=0.01$. When the regime of convergence is reached, the parameter $\epsilon _{n}^{{\rm max} }$ is decreasing rapidly and the ratio ${{q}_{n}}=\epsilon _{{\rm max} }^{(n+1)}/\epsilon _{{\rm max} }^{(n)}$ is converging to a constant value.

	${{J}_{2}}/J=0.5$		${{J}_{2}}/J=0.01$
${{n}_{B}}=10\times {{2}^{n-1}}$	$\epsilon _{{\rm max} }^{(n)}$	qn	$\epsilon _{{\rm max} }^{(n)}$	qn
n = 2	0.197	0.419	0.693	0.847
n = 3	0.087	0.471	0.588	0.713
n = 4	0.041	0.496	0.419	0.516
n = 5	0.020	0.496	0.216	0.472
n = 6	0.010	0.497	0.102	0.489
n = 7	0.005	0.499	0.049	0.491

Next, the numerical factor of the error is obtained using equation (F.9). Figure F.2 displays the Berry conductivity for three successive iterations n_B = 80, 160 and 320 with error bounds. Here the parameter n_R = 200 is chosen such that the statistical error is negligible. Otherwise, one should sum the statistical error to the error due to the grid. The error bounds are becoming smaller at each iteration by a factor 2, and the conductivity converges rapidly with an error already of $0.03\ {{\mathfrak{C}}^{(n}})$ for n_B = 160 and $0.015\ {{\mathfrak{C}}^{(n)}}$ at n_B = 320.

We have finally applied the algorithm on a more challenging case of the Haldane model with the parameters ${{J}_{2}}/J=0.01,\beta /J=0$. In this parameter regime, the band structure exhibits only a small energy gap between the two bands, and the Berry curvature is very peaked at the Dirac cones. However, the FHS algorithm is already working at n_B = 10. As can be inferred from the right column of table 1, the regime of convergence is only reached for a finer choice of the grid than in the case ${{J}_{2}}/J=0.50$, with a slower decrease of $\epsilon _{{\rm max} }^{n}$ at the beginning. However, at n = 4 the system enters a convergence regime; and also in this case, one can associate with the Berry conductivity upper error bounds due to the fine discretization of the Brillouin zone.