Efficient algorithm to compute the Berry conductivity

We propose and construct a numerical algorithm to calculate the Berry conductivity in topological band insulators. The method is applicable to cold atom systems as well as solid state setups, both for the insulating case where the Fermi energy lies in the gap between two bulk bands as well as in the metallic regime and interpolates smoothly between both regimes. The algorithm is gauge-invariant by construction, efficient and yields the Berry conductivity with known and controllable statistical error bars. We apply the algorithm to several paradigmatic models in the field of topological insulators, including Haldane's model on the honeycomb lattice, the multi-band Hofstadter model and the BHZ model, which describes the 2D spin Hall effect observed in CdTe/HgTe/CdTe quantum well heterostructures.


I. INTRODUCTION
Topological insulators (TI) are a topological state of quantum matter which constitutes a new paradigm in condensed matter physics [1][2][3][4] . These recently discovered new materials exhibit unique fascinating properties such as currentcarrying surface or edge states that are strongly protected against perturbations in either the bulk or the surface of the material 5-10 and non-standard exchange statistics of quasiparticle excitations, which offer potential applications in the context of quantum computation [11][12][13] .
The question what happens in topological insulators when the Fermi energy does no longer lie inside the gap between two energy bands, is by no means rhetoric 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 2 Se 3 and Bi 2 Te 3 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 one 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 2 Se 3 and Bi 2 Te 3 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 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 2e 2 /h.
Complementary to solid-state realizations, cold atoms in optical lattices have been proposed as a realistic plat-form to experimentally explore the new physics of TIs under controllable conditions [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] . In particular, in these systems the Fermi energy can be controlled directly by the filling of atoms in the lattice. 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 plays as well 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 ρ x y 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 the complications arising from the fact that there are three main mechanisms that influence the electronic motion and can give rise to an AHE. 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 scattering that affects other AHE mechanisms.
Understanding of 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 FIG. 1. a. Generic energy spectrum of a system with an energy gap ∆. In the displayed situation the Fermi energy falls into the first energy band and defines the Fermi surface as the equipotential energy line at E α (k) = E F (solid line). The projection of the energy dispersion of the first band is shown as a color-coded plot in the horizontal k x − k y -plane. b. For the numerical calculation of the Berry conductivity, the Brillouin zone is discretized by a finite grid. Momentum space plaquettes with energies E α (k) entirely below (above) the Fermi energy contribute entirely (not at all) to the Berry conductivity, whereas plaquettes which cut the Fermi surface contribute partially. c. Schematic summary of the numerical algorithm to calculate the Berry conductivity: After fixing the discretization grid of momentum space and calculating the Berry curvature contributions by means of the FHS algorithm for each plaquette of the Brillouin zone, a classical Monte Carlo sampling method is used to determine the weights with which the individual plaquettes contribute to the conductivity. Statistical uncertainties in the sampling process result in controlled and statistical errors in the Berry conductivity. 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 and 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, which 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 to the Chern number according to the 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 Fig. 1), discuss its ingredients and show that it is gauge invariant and efficient (Sec. II).
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 since when the Berry conductivity is not integervalued, errors due to approximations need to be under control in order to distinguish two different values of the conductivity observable, so that one 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 (Sec. III A). 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 (Sec. III B). 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 BHZ model 16 which is a realistic model which captures the physics of 2D spin Hall effect present in systems such as the above-mentioned CdTe/HgTe/CdTe quantum well compounds (Sec. III C). We conclude with a short summary and a discussion of possible future extensions of the presented method (Sec. IV).

A. 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 both for the insulating case where the value of the Fermi energy lies in the gap between two bands, as well as the generalized result for the situation in which the Fermi energy falls 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, The Chern numbers C α are integer-valued topological invariants, defined in terms of the integral of the Berry curvature F α x y (k) over the whole Brillouin Zone (B.Z.) 39,42 : The latter is expressed by the exterior derivative of the Berry connection where u α (k) is the eigenvector corresponding to the energy band E α (k).
In the case that the Fermi energy does not lie in an energy gap between bands, as schematically shown in Fig. 1a, the intrinsic Hall conductivity generalizes to 36,43 with where Θ(E ) denotes the Heaviside function. Thus, the conductivity is the sum of the integer-valued Chern numbers corresponding to fully-occupied energy bands below the Fermi energy and a non-quantized contribution which 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, as e.g. in two-band systems, the expressions for the eigenenergies and eigenvectors of the bands are given in explicit form, and hence the Chern values C α can be calculated analytically. In general, however, the system Hamiltonian cannot be diagonalized analytically and an efficient numerical method to compute the Chern values is needed.

B. Construction and properties of the algorithm
The algorithm we propose to numerically compute the Chern values of Eq. (12) and thereby the Berry conductivity of Eq. (4) is based on a series of controlled approximations: First, we discretize the two-dimensional Brillouin zone by a finite n B × n B grid of small plaquettes at discrete momenta k l (see Fig. 1b and Appendix A for details), so that the integral over the (partially filled) band becomes with the Berry curvature contribution from a small two-dimensional plaquette of size ∆ k x ∆ k y , and the weighting factors The weights p α l (E F ) correspond to the partial area of the plaquette, which is covered by the Fermi sea, thus p α l (E F ) = 0 (p α l (E F ) = 1) for squares with energies completely above (below) the Fermi energy E F , and 0 < p α l (E F ) < 1 for momentum space plaquettes which are cut by the Fermi surface (see Fig. 1b). The choice 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.
The key of the numerical algorithm is now to evaluate reliably and under controlled approximations the discretized sum of Eq. (6), whose value converges to Eq. (12) for increasingly finer grids.
(i) Gauge-invariant calculation of the Berry curvature: To numerically calculate the Berry curvature contributions F α x y,l we employ a numerical algorithm proposed by Fukui, Hatsugai and Suzuki 44 (FHS algorithm). It is highly efficient and the discrete sum 1/(2πi ) {k l }F α x y,l converges rapidly to the correct integer-valued Chern numbers C α , 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 strengthF α x y,l calculated with the FHS method.
(ii) Efficient estimation of the weights p α l (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 k l , we generate n R uniformly distributed random points k R and compute E α (k R ) which takes the discrete values zero or one.
Based on the latter we define the estimator for the weighting factors p α l (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 that the value of the estimatorsp α l (E F ) is bounded between zero and one. However, it is clear that the statistical error will be largest for partially contributing plaquettes withp α l (E F ) ∼ 1/2, whereas the uncertainty inp α l (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 inp α l (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 which depend on the actual value of the Fermi energy. The quantityp α l (E F ) is the estimator of the fixed though unknown parameter p of a binomial distribution 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 < 1 this allows one to derive a confidence interval [p α l ,min , p α l ,max ] for p α l (E F ), called the Wilson interval 47 with modified boundary conditions. This means that with a probability 1 − the "true" value p α l lies in this interval. The key point is that the width of this interval depends on the actual value of the estimatorp α l and is typically significantly smaller than the trivial upper bound of one. After symmetrizing the interval by taking the maximum ∆p α l (E F ) = max(p α l − p min , p max −p α l ), each momentum space plaquette of the grid is associated to a probability valuep l (E F ) ± ∆p l (E F ) with a confidence of at least 1 − .
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 ofp α l directly translates into an uncertainty in the Berry conductivity contributions, Finally, the estimated Berry conductivity is given bỹ with an error ±∆C α (E F ) of with confidence of at least 1 − . We remark that controllable error bars in combination are particularly important and valuable outside of the insulating regime, i.e. where the Fermi energy cuts a partially filled energy band, as in this case the Berry conductivity is not quantized and can assume continuous non-integer values.

III. PRACTICAL APPLICATION OF THE ALGORITHM
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 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.

A. The Haldane model
The model proposed by Haldane in 40 is a tight-binding Hamiltonian of spinless fermions on a honeycomb lattice, with dynamics governed by nearest-neighbor (N.N.) realvalued hopping term of amplitude J and an imaginary nextto-nearest neighbor (N.N.N.) hopping term J 2 (see Fig. 2a). In addition, the fermions are exposed to an onsite staggering potential β, which induces a chemical potential difference between nearest-neighbors 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 .
The Hamiltonian of the system is given by Here, c † i and c i are fermionic creation and destruction operators, ν i j = sgn[(d 1 × d 2 ) z ] and s φ,ψ = ±1. The vectors d 1 and d 2 are oriented along the bonds of the hexagonal unit cell, as shown in Fig. 2a. The model can be readily solved by rewriting the Hamiltonian in terms of two-site basis cells (φ, ψ) (see e.g. 48 ) such that the hexagonal lattice becomes a triangular lattice of (φ, ψ) cells. In the Fourier space the Hamiltonian is then given by 40 is expressed in terms of the vectors between nearest neighbors δ 1 , δ 2 and δ 3 and f (k) = sin[a 1 ·k]+sin[a 3 · k]+sin[(a 1 +a 2 )·k] is expressed in terms of the lattice vectors a 1 and a 2 as shown in Fig. 2 and defined in Appendix C.
Diagonalization of the Hamiltonian readily yields the two-band energy spectrum which is shown in Fig. 2b. For β = J 2 = 0, the Hamiltonian corresponds to pure nearest-neighbor hopping of fermions with the characteristic spectrum exhibiting the two inequivalent Dirac cones 49,50 . A non-zero staggering potential β = 0 induces an imbalance of the fermion density on φ and ψ lattice site. The formation of a charge-density-wave phase is associated to the opening of a topologically trivial insulating gap in the spectrum. On the other hand, a strong enough N.N.N. hopping term J 2 opens a topologically non-trivial energy gap that signals the transition of the system into a AHE phase characterized by a non-zero Chern number. The size of the energy gap is determined by the formula ∆ = 2|β − 3 3J 2 |, and for |β| < 3 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 Fig. 1c -to the Haldane model. To this end, we start by fixing the Hamiltonian parameters to J 2 = 0.1J , β = 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 strengthF x y for each plaquette (step 2); the result is shown in the right column of Fig. 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 α (k R ) (step 4). Once the Fermi energy is fixed (step 5), here to a value of E F = −1.5J values of n B , illustrating how an increasingly finer grid of the Brillouin zone leads to an increased resolution and numerical precision.
Finally, the estimated weightsp α l (E F ) and the Berry curvature contributionsF x y,l are combined to calculate the Berry conductivity (step 8) according to Eqs. (11) and (12) with an associated error bar (step 9) as given by Eq. (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 Fig. 4a: starting from low conductivity values at the bottom of 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.
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 Fig. 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 few plaquettes -a fact that the algorithm signals in the form of larger error bars of tem 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 σ B e (E F ) = (e 2 /h) 3 3J 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 Fig. 5.

B. 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 Φ. Only very recently, several groups have achieved to observe the characteristic physics including the fractal spectrum known as Hofstadter's butterfly in graphene superlattice systems [54][55][56] . This is complementary to ongoing experimental efforts to realize theoretical ideas 57 on how to implement the fermionic Hofstadter Hamiltonian with cold atoms in optical lattices [58][59][60][61] .
The Hamiltonian in second-quantized form is given by where the sum is over nearest neighbor sites (see Fig. 6) and the phase factor exp(i θ i j ) 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 Φ = 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. 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. Fig. 7 shows the numerical results for the Berry conductivity for different values of the flux per plaquette (Φ = 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 Φ = 1/3 the Berry conductivity interpolates monotonically between the gap plateau values, for Φ = 1/5 the conductivity displays an interesting feature for Fermi energy values in the second band: instead of showing a 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 pronounced, in the third band of the spectrum for Φ = 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.

C. The BHZ model
In 2005, it was suggested that the quantum spin Hall effect (QSHE) could possibly be observed in graphene 51,62 , which 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 Z 2 topological invariant. The physics can be described by an effective Hamiltonian valid close to the Γ point, derived by Bernevig, Hughes and Zhang (BHZ model) 16,63 . The Hamiltonian is a given by 4 × 4 matrix in momentum space, , where 1 is the two-dimensional identity matrix, σ i denote the Pauli matrices and 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 and it makes thus 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 8a shows the energy spectrum that exhibits a small gap of 0.01eV , which renders the computation of the Berry conductivity in the non-insulating regime more demanding. Figure 8b -d shows the numerical results for the Berry conductivity for increasingly finer grids of the Brillouin zone.
Whereas even for the roughest grid studied (n B = 40) the algorithm correctly captures the qualitative behavior and the conductivity minimum value value of -1 for the Fermi energy lying in the shallow energy gap. However, as signaled by considerably large error bars, only few plaquettes contribute large values of Berry curvature to the conductivity. Thus, finer grids (see Fig. 8c and d with n B = 160 and n B = 320) are required to quantitatively correctly describe the conductivity behavior in the vicinity of the gap. 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 coarsegrained grids, this can be helpful to predict observations in the case of restricted experimental resolution, e.g. originating from limited sizes of optical lattices for cold atoms, or finite temperature constraints in solid state experiments.
Finally, we remark that the BHZ model is an effective model valid close to the Γ 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.

IV. CONCLUSIONS AND OUTLOOK
In this work we have proposed and constructed a numerical algorithm to calculate the Berry conductivity in topological band insulators. The algorithm works for both the insulating case where the Fermi energy lies in the gap between two bulk bands as well as the situation where 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 which 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][69][70] . Recent experiments in which topological 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 and can predict some interesting features within the energy band. In fact, there have been proposed several ways to measure characteristic signatures of topological quantum phases in systems of cold atoms [71][72][73][74][75] .
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, it is an interesting question 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][77][78][79][80] . The continuous Brillouin Zone is discretized by a twodimensional 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 with 0 ≤ i , j ≤ n B − 1, The lattice field strengthF α x y (k l ) of band α on the grid is then defined in terms of the link variable U µ (k) as If the admissibility condition |F α x y (k l )| < π is satisfied 44,46 , the lattice gauge theory corresponds to the continuous gauge theory 44,46 and one can write: Based on these Berry curvature contributions, the Chern number can be computed as

Appendix B: Choice and the computation of the statistical error
In this section, we present the concept and the details of a confidence interval (C.I.) to characterize the statistical uncertainty of the estimated weightsp α l , as defined in Eq. (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 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 The goal is to associate a C.I. of a width much smaller than one to the estimated valuep, such that the true value p lies with a probability 1− 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 − . As a guiding principle, a "good" C.I. is an interval with p cov 1 − . On the contrary, for p cov < 1 − , the C.I. is "bad" as the statistical "guaranteeing functionality" of the interval fails. The other case p cov > 1 − 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 − .
The construction of 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 N , in our case: The central limit theorem and the definition of the C.I. of a normal distribution with an expected probability 1− permits us to write the C.I. ofp l as a self-consistent equation in terms of p: where z 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 for p = 1/2. This relation highlights the typical 1/ 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 does not longer depend on the estimated valuep l itself, it does not satisfy our requirement. It will have a coverage probability p cov > 1 − 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) p l (1−p l ) in Eq. (B3), thereby replacing the unknown "true" value by the estimator value, so that This C.I. is known as the Wald interval 47 . Despite its simplicity, this convention suffers from several problems: for p l 0 orp l 1, the Wald interval shrinks to zero, implying a bad a 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 Fig. 9a 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 10000 samples. One notices at first glance the tendency of the curve to lie below the expected value of 1 − . 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 a 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 cov to a poor p cov value as it is the case for instance around p = 0.8 in the shown example. The couple (p, n R ) defines the lucky/unlucky numbers. In Fig. 10a we fix the value p = 0.25 and vary the value of n R . Here one observes also significant fluctuations that are only stabilized at larger values of n R . This effect becomes is even more striking at small p, as illustrated in the Fig. 10 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 cov until reaching n R = 423 where p cov suddenly drops from 0.94 to 0.78. We thus exclude the Wald interval as a candidate to construct the C.I. for thep α l estimators in our algorithm.
Most of the mentioned problems can be avoided if the approximation p(1 − p) p l (1 −p l ) is not applied in Eq. (B3). Instead, one can exactly solve Eq. (B3), which is a quadratic equation forp. This yields the so-called Wilson interval 47,82 As illustrated in the Fig. 9b, the Wilson interval is much more stable and the coverage probability is oscillating around the value 1− . Figure 10 b. shows that the Wilson interval reaches rapidly and in a stable way the expected value 1 − . The only problem still to be cured is at the boundaries, at p-values around zero or one, where the coverage probability drops. Figure 10 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 ofp. 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 ef- fect: we use the following replacement: 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. aroundp by choosing a width which corresponds to twice the value of max{p max −p,p − p 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.
In this section, we illustrate the importance of an appropriate momentum space resolution, parametrized by the discretization number n B . Figure 12 presents a zoom of Fig. 4c 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.33. There, the behavior of the estimator of the Berry curvature becomes jerky. This is a characteristics which shows up when some few momentum-space plaquettes have an important Berry curvature contribution. The error bars are signal this effect. The situation improves for increasing values of n B : the curves converging to one sharp curve, showing that the main contribution of the Berry curvature is stems from states with an energy close to zero, and the error bars decrease significantly. FIG. 13. Numerically estimated Berry conductivityσ Be with error bars for two values of the number of random points n R = 20 and 160 for a system with J 2 = 0.1J , β = 0.5J and n B = 20. As expected, the computation with n R = 160 is much more precise, resulting in significantly smaller error bars. Note that as desired the n R = 160 curve is entirely comprised in the region spanned by the error bars of the computation with n R = 20.
Another way to reduce size of the error bars is to increase the value of n R , the number of random points used to computep l in each plaquette. Figure 13 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.

Appendix E: Importance of the choice of the error bars
In this section, we compare the Wilson interval with modified boundaries with the C.I. defined in Eq. (B4) by examining the final error interval obtained in the Haldane model using both methods. We work here with J 2 = 0.1J , β = 0.5J such that the Berry curvature is really sharp and localized. The Figure 14 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 Fig. b., is constant and independent of the value of the Fermi en-