Ab initio analysis of the topological phase diagram of the Haldane model

We present an ab initio analysis of a continuous Hamiltonian that maps into the celebrated Haldane model. The tunnelling coefficients of the tight-binding model are computed by means of two independent methods - one based on the maximally localized Wannier functions, the other through analytic expressions in terms of gauge-invariant properties of the spectrum - that provide a remarkable agreement and allow to accurately reproduce the exact spectrum of the continuous Hamiltonian. By combining these results with the numerical calculation of the Chern number, we are able to draw the phase diagram in terms of the physical parameters of the microscopic model. Remarkably, we find that only a small fraction of the original phase diagram of the Haldane model can be accessed, and that the topological insulator phase is suppressed in the deep tight-binding regime.


I. INTRODUCTION
The Haldane model [1] is a celebrated lattice model describing a Chern insulator [2], characterized by the presence of quantum Hall effect (QHE) [3] in the absence of a macroscopic magnetic field. Conceptually, the Haldane model stands at the heart of the tremendous advances in the field of topological condensed matter physics, as the mechanism for a non-trivial band topology presented by Haldane is realized in actual materials via the intrinsic spin-orbit interaction of topological insulators [4,5]. These concepts are also relevant for the physics of ultracold atoms in optical lattices, as these systems represent a powerful platform for simulating solid-state physics [6]. Mott insulators [7,8], bosonic superfluids [9] or graphenelike honeycomb lattices [10][11][12][13][14][15][16] are among the many systems that have been emulated by this technique. Interestingly, an effective experimental realization of the Haldane model has been recently reported in Ref. [17].
In his original work, Haldane constructed a discrete tight-binding model for a non-centrosymmetric honeycomb lattice in the presence of a vector potential A(r), with vanishing total flux through the unit cell. The key feature of the model is that, even in absence of a macroscopic magnetic field, the time-reversal symmetry is broken due to the presence of the gauge field A(r). This, in turn, implies that the next-to-nearest neighbour tunnelling coefficient t 1 becomes a complex number. Haldane showed that the properties of the system depend on the interplay between the phase acquired by t 1 and the effect of parity breaking, affecting the topological phase diagram of the model [1].
Considering the above, the knowledge of the dependence of the phase acquired by t 1 on the applied vector potential field becomes a crucial element for drawing the topological phase diagram. For this purpose, it is common practice [1,18] to make use of the so-called Peierls substitution, whereby the effect of A(r) is effectively included by the replacement t 1 → t 1 exp (i(e/h) A(r)dr) [19]. However, in a recent work [20] we showed that the Peierls substitution is actually wrong whenever the vector field A(r) has the same periodicity of the underlying lattice, as it is the case of the Haldane model by construction. In that work, we analyzed the parity invariant case by presenting two independent approaches for calculating the tight-binding parameters of the model: one based on the maximally localized Wannier functions (MLWFs), the other on a closed set of analytical expressions in terms of the energy spectrum at selected high symmetry points in the Brillouin zone (BZ).
In the present work, we extend the previous analysis to the general case in which both inversion and timereversal symmetry can be broken. We show that the two approaches considered provide a remarkable agreement even in the presence of parity breaking, allowing for a precise determination of the tight-binding parameters of the model. By combining these results with the numerical calculation of the Chern number, we are able to redraw the topological phase diagram of the Haldane model in terms of the physical parameters of the microscopic model. Interestingly, we find that only a small fraction of the original phase diagram can be accessed, and that the topological insulator phase shrinks dramatically as the system becomes more and more tight-binding. In addition, we find that the gap closing at the topological phase transition does not take place exactly at one of high symmetry points of the BZ, but in a close-by point. The reason is that the complex tunneling between ho-mologous sites are no longer degenerate in the presence of parity breaking, contrarily to what it is assumed in the Haldane model.
The paper is organized as follows. In section II we introduce the microscopic continuous Hamiltonian used in this work and review the formal steps needed to derive the corresponding tight-binding model. Some general properties of the Haldane model are also recalled. Then, in section III we present the two approaches employed for calculating the tight-binding parameters, and discuss how the breaking of time-reversal and/or parity affects their behavior. In section IV we analyze the topological phase diagram, both in terms of the parameters of the tight-binding model and of the physical ones. Concluding remarks are drawn in section V. Finally, in the appendices we present an analysis of the spread functional of the MLWFs (Appendix A) and additional remarks on the numerics (Appendix B).

II. SETUP OF THE HALDANE MODEL
In this section we present a systematic derivation of the Haldane model starting from the continuous Hamiltonian proposed by Shao et al. [18] in the context of cold atoms trapped in optical lattices (see also [20]). The method discussed here is general and suited to map a generic continuous Hamiltonian to its corresponding tight-binding model [21,22].

A. The continuous Hamiltonian
The general form of a continuous Hamiltonian in the presence of a scalar lattice potential V L (r) and a vector potential A(r) is with r = (x, y) in case of a two dimensional system, as we shall consider here. The potential V L (r) is chosen in order to generate a periodic structure with two minima per unit cell, forming a honeycomb lattice [10,18]: Above, E R =h 2 k 2 L /2m is the recoil energy, k L denotes the laser wavevector, s is a dimensionless parameter representing the strength of the potential in units of 3e y ) are the basis vectors in the reciprocal k space, and χ A is a parameter related to the breaking of the parity symmetry. In particular, χ A = 0 corresponds to the inversion symmetric case, where the two minima in the unit cell are degenerate. On the other hand, for χ A = 0 parity is broken and the minima are no longer degenerate. The unit cell (shown in Fig. 1) is generated by the direct lattice vectors a 1,2 = (2π/3k L )(e x ∓ √ 3e y ). We also define the lattice vector a 3 = a 1 + a 2 , that will be useful later on.
We now turn to the vector potential contained in the Hamiltonian (1). As already mentioned, we employ the form proposed by Shao et al [18], namely with ∇ · A(r) = 0 (Coulomb gauge). The flux of the corresponding magnetic field B = ∇×A across the unit cell is null [18], as required for the Haldane model. In the following analysis, the only variable parameter entering the expression for the vector potential A(r) is its amplitude α. Notice that for α = 0 the system is symmetric under time-reversal, whereas this is not the case for α = 0.

B. The tight-binding model
The continuous Hamiltonian (1) can be mapped onto the tight-binding Haldane model [1,18] by following the general procedure discussed in [20][21][22]. The starting point is the (single particle) many-body Hamiltonian defined byĤ withψ a field operator for bosonic or fermionic particles. Then, when the wells of the lattice potential are deep enough, the field operator can be conveniently expanded in terms of a set of functions w jν (r) localized at each minimum:ψ (r) ≡ jνâ jν w jν (r).
Above, ν is a band index andâ † jν (â jν ) is the creation (destruction) operator of a single particle in the jth cell, satisfying the usual commutation (or anti commutation) rules following from those of the fieldψ.
As already mentioned in the introduction, we shall use the maximally localized Wannier functions (MLWFs) for composite energy bands [23] as basis functions. The MLWFs are defined as linear combinations of the Bloch eigenstates ψ ν k (r), where S B represents the volume of the first BZ, and U ∈ U (N ) is a gauge transformation that is obtained from the minimisation of the Marzari-Vanderbilt spread functional, Ω = ν r 2 ν − r 2 ν [23]. We remark that the presence of the vector potential may significantly affect both the Bloch eigenfunctions ψ ν k (r) and the unitary matrices U νν (k) [20]. A thorough analysis of the MLWFs is given in section III A.
In the following, we shall consider the contribution of the first two Bloch bands only, namely ν, ν = 1, 2. This is sufficient for constructing the lowest lying ML-WFs localized at the two lattice sites A and B inside the unit cell. Then, by considering the following transformation from coordinate to reciprocal space, b νk = (1/ √ S B ) j e −ik·R jâ jν , the reduced tightbinding Hamiltonian can be written aŝ where the 2 × 2 matrix h νν (k) is given by with R j = {j 1 a 1 + j 2 a 2 j 1 , j 2 = 0, ±1, ±2 . . . } and j labels the unit cell. We remark that the eigenvalues of h νν (k) coincide with the exact Bloch energies n (k) (n = 1, 2) if the full expansion of neighbouring coefficients R j is retained. When the system is in the tightbinding regime (s > ∼ 5) [21], it is convenient to truncate the series by retaining only a finite number of matrix elements w 0ν |Ĥ 0 |w jν , with the eigenvalues of h νν (k) still being a good approximation of the exact energies. We note that since the functions w jν (r) are in general complex (see section III A), we may expect the matrix elements to be complex as well.
By truncating the tight-binding expansion in Eq. (8) to the next-to-nearest neighbours (see Fig. 1), the matrix h νν (k) can be written as the sum of three terms: The first term corresponds to the on-site energies which are real quantities by definition. The next term, h νν , contains only off-diagonal elements corresponding to the hopping between the three nearest-neighbour sites. Although the basis functions w jν (r) are complex, these three tunnelling amplitudes can be chosen to be real by means of a suitable global gauge fixing, as they are all equal thanks to the symmetries of the system (see Fig.  1). Taking this into consideration, and defining The rotational symmetry implies that next-to-nearest tunnelling amplitudes t1 along the same direction are conjugate pairs (solid and dashed lines in red); from the latter follows the equivalence of the hopping amplitudes separated by 2π/3 radians. When sites A and B are degenerate, the system is also invariant under rotations by π radians around the centre of any elementary cell.
we can write Finally, the term h (2) νν (k) corresponds to the six next-tonearest neighbours hopping between homologous sites. By taking into account all the symmetries of the full Hamiltonian of Eq. (1), these tunnelling coefficients can be compactly written as The above equation explicitly shows that t ± 1ν contains two distinct complex phases ±ϕ ν for each site type (ν = A, B). Then, using Eq. (13) and after some algebra, we can write The above expressions allow to cast the matrix h νν (k) in the following compact form, where we have defined By expanding the Hamiltonian on the basis of the 2×2 identity matrix, I, and of the Pauli matrices, σ i , Eq. (15) can be rewritten as [18] where the coefficients h i (k) are given by the following expressions: with the vectors s i being those joining the three nearest neighbours of type A(B) to a given site of type B(A) [18].
In the last expression we have also fixed, without loss of generality, E A = and E B = − . Then, the tight-binding energies are readily found as

The Haldane model and the Peierls substitution.
At this point, further approximations are required in order to recover the original model proposed by Haldane [1], namely |t 1A | = |t 1B | ≡ |t 1 | and ϕ A = −ϕ B ≡ ϕ. We shall refer to this configuration as the "simplified parameter setup" (SPS). We note that the corresponding model contains only four parameters, namely , t 0 , |t 1 | and ϕ. In section III C, we shall provide numerical evidence showing that in the tight-binding regime the difference between |t 1A | and |t 1B | is negligible, thus justifying the SPS. The second condition is not strictly verified (again, see section III C), but one can consider a sort of effective model by defining a single phase, ϕ ≡ (ϕ A − ϕ B )/2. Therefore, in the SPS the terms h 0 and h 3 of Eqs. (18) and (21) simplify to The above equations, together with Eqs. (19) and (20), correspond to the definition of the original Haldane model, see [1,18]. As already mentioned in the introduction, the original model proposed by Haldane is constructed by means of the so-called Peierls substitution [1,18]. This consists in assuming that the complex phase of the tunnelling coefficient t ij is given by the integral of the vector potential along the straight path joining sites i and j, i.e. t ij → t ij exp{i(e/h) j i Adr}. In the present case, the Peierls prediction for the complex phase would be [18] This value will be used later on, in Secs. III and IV, for comparison with the results of the two approaches discussed in this paper. Here, we can anticipate -as thoroughly discussed in [20] -that the above result is definitely incorrect, owing to the fact that the Peierls substitution is justified only when the vector field A(r) is slowly varying on the scale of the lattice [24]. In fact, this condition is explicitly violated in the Haldane model, where A(r) has the same periodicity of the lattice.

C. General features of the Haldane model
Before proceeding, let us recall some general features of the Haldane model [1], corresponding to the SPS approximation. This model is characterized by the presence of Dirac points located at the vertices k D of the first BZ, where the dispersion law takes the relativistic form (k) = √ m 2 c 4 + c 2 k 2 . They can be divided into two inequivalent classes, corresponding e.g. to k ± D = ±(1, 0)k L (often referred in the literature as K and K ; in the following, we will always refer to these two as inequivalent Dirac points, for simplicity). In the presence of time-reversal and inversion symmetry (namely, for α = 0, χ A = 0), the two energy bands become degenerate at the Dirac points, whose position is given by the solution of z(k D ) = 0 or h 1 (k D ) = h 2 (k D ) = 0 (this holds at any order of the tight-binding expansion). In fact, in this case both and ϕ vanish, yielding h 3 (k) = 0, f A = f B , and f − = 0, so that the energies ± (k D ) are degenerate.
When both time-reversal and inversion symmetry are broken (α = 0, χ A = 0), two inequivalent energy gaps form at the Dirac points: The closure of one of them indicates a topological phase transition, where This equation identifies the well known boundary between the normal and topological insulator phases with Chern numbers C = 0 and C = ±1 in the Haldane model [1]. We remark that in the general model the gap closing does not take place exactly at k ± D , but in a close-by point. This is due to the fact that, when breaking of parity is included self consistently, the tunnelling parameters t 1A and t 1B are no longer degenerate, contrarily to what it is assumed in the Haldane model (cf. the SPS approximation).

III. CALCULATION OF THE TIGHT-BINDING PARAMETERS
In this section we discuss two independent methods for calculating the tight-binding parameters for arbitrary values of the physical parameters s, α and χ A . The first method is based on the ab initio calculation of the maximally localized Wannier functions (MLWFs) [23,25], which we already employed in [20][21][22] in different lattice geometries. This approach gives direct access to the whole set of parameters , t 0 , |t 1A |, |t 1B |, ϕ A and ϕ B . The second approach relies instead on analytical expressions in terms of the energy spectrum, as discussed in [20]. The latter is here extended to the case of parity breaking.
We remark that the approach based on the MLWFs corresponds to the ab initio definition of the parameters, and is therefore model-independent. Instead, the second method depends on the specific form of the tight-binding Hamiltonian. However, it does not require the calculation of any set of Wannier functions since only the spectrum of the continuous Hamiltonian is needed. Notably, the two methods present a remarkable agreement in the whole range of parameters considered here.

A. Maximally localized Wannier functions
The MLWFs, which have been defined in Eq. (6), represent a powerful tool that is largely employed in condensed matter physics [25]. By construction, the ML-WFs are the basis functions with the maximal degree of localization in coordinate space, allowing to construct tight-binding models that accurately reproduce the properties of the continuous Hamiltonian, with a minimal set of tunnelling coefficients. In addition, the MLWFs permit a very fine sampling of the reciprocal space thanks to the so-called Wannier interpolation technique [25]. This point is very important for our purposes in this work, as the determination of the Chern number requires a high density of points in k-space [26,27].
The MLWFs are computed by means of the standard implementation of the WANNIER90 package [28] (see also Appendix B). The resulting functions are complexvalued when α = 0. This feature is in agreement with the analysis of [29], where it was shown that, in general, MLWFs cannot be constructed as real functions when the time-reversal symmetry is broken (see also [30]). In the context of the Haldane model, the imaginary part of the MLWFs plays an essential role, since it determines the complex phase acquired by the next-to-nearest tunnelling coefficient. In turn, the complex phase directly affects physically meaningful quantities, such as the spectrum or the topological phase diagram [1].
In Fig. 2 we illustrate an example of the real-space structure (note the logarithmic scale) of the real and imaginary parts of a MLWF for sublattice ν = A, located at the origin j = 0, for s = 5, α = 0.1 and χ A = 0. The structure of the real part is very similar to the one in the pure honeycomb lattice [21], namely it is highly localized around the origin, with appreciable contribution around the neighbouring lattice sites. In average, the imaginary part is two-three orders of magnitude smaller than the real part. It is particularly interesting to observe that the imaginary part is null at the interstitial region between nearest neighbours, while it becomes maximum along the path joining next-to-nearest neighbours. These properties hold in the whole range of parameters considered in this work.
An analysis of the spread functional of the MLWFs as a function of the vector potential amplitude has been included in Appendix A.

B. Analytical expressions from the spectrum
In this section we derive a closed set of analytical expressions in terms of the energy spectrum at selected high symmetry points in the BZ. This is done in the framework of the SPS discussed in section II B 1, corresponding to the standard formulation of the Haldane model [1,18]. As we shall see below, the approximations of the SPS are well justified in the tight-binding regime. The model is therefore given in terms of four parameters, namely , ϕ, t 0 and |t 1 |. We remind that measures the difference between the on-site energies E A and E B , and it is therefore associated to the breaking of parity, whereas the breaking of the time-reversal symmetry corresponds to ϕ different from zero. It is also worth recalling that the parameters of the underlying continuous Hamiltonian that control the breaking of parity and time-reversal symmetry are χ A and α, respectively. In particular, χ A = 0 gives = 0 whereas α = 0 implies ϕ = 0. We begin by noting the following relations at k = 0: Similarly, at the Dirac points k ± D we have z(k ± D ) = 0.
Next, let us define the bandwidths Recalling the expression for the gap at the Dirac points in Eq. (26), one can easily derive the following relations: Due to the symmetries of the system, we can consider ≥ 0 and ϕ ≥ 0 without loss of generality. Focusing first on the region with > 3 √ 3|t 1 | sin ϕ (corresponding to the normal insulator phase), after some algebra one finds the following set of formulas: Similarly, in the region with < 3 √ 3|t 1 | sin ϕ (corresponding to the topological insulator phase), we find the following expressions: The solutions in a generic case with < 0 or ϕ < 0 can be obtained from symmetry considerations, by exchanging the role of the two basis points A, B and/or of the two inequivalent Dirac points k ± D .

C. Numerical results
In this section we present a comparison of the two methods described in Secs. III A and III B for the calculation of the tight-binding parameters. In addition, we also analyze the accuracy of the assumptions of the SPS (Sec. III B) based on the tunneling coefficients extracted from the MLWFs.
Let us begin by analyzing Fig. 3, where we compare the tunnelling coefficients calculated from the MLWFs with those calculated from the analytical formulas of Eqs. (37)-(40), valid for the normal insulator regime, and Eqs. (41)-(44), valid for the topological insulator regime, which is depicted by the grey shaded area in the figure.
Results are shown as a function of α for fixed values s = 5 and χ A = 0.001, since essential features are unaffected by s and χ A . In the case of the MLWFs, we have plotted the averages ϕ = (ϕ A − ϕ B )/2 and |t 1 | = (|t 1A | + |t 1B |)/2 in order to allow comparison with the analytical formulas, which have been derived in the context of the SPS (see section III B) Overall, Fig. 3 shows a very good agreement between the two methods for all the tunnelling coefficients, in all regimes. Furthermore, it is interesting to note that the two different solutions represented by the set of Eqs.  To put in other words, the solution of one set of equations on one side represents a smooth continuation of the solution of the other set of equations in the other side, and viceversa. Provided that one chooses the right solution, the calculated values agree very well with those of the MLWFs, as already said. In addition, Fig. 3(d) reveals an extremely important feature that was absent in the original Haldane model: the phase ϕ is limited by a maximal value. This behavior, that was already found in the parity-symmetric case [20], implies that ϕ can only access a restricted range of values, therefore limiting the physically accessible region of the phase diagram. This feature will be crucial for the analysis presented in the next section, where we shall redraw the topological phase diagram in terms of the physical parameters -α, χ A and s -of the underlying continuous Hamiltonian. Next, we proceed to test the accuracy of the assumptions of the SPS approximation (Sec. III B) based on the tunneling coefficients calculated from the MLWFs. This is done in Fig. 4, where we compare the relative deviations from the average values of the phase, 1 − ϕ A,B /ϕ, and of the magnitude of the next-to-nearest tunneling coefficient, 1 − |t 1A,B |/|t 1 |, for χ A = 0.001, s = 5. This figure demonstrates that the maximum relative deviation in both cases is below ∼ 1%. We have verified that this holds for all values of s and χ A considered here, thus justifying the assumptions of the SPS approximation in the whole range of parameters. Apart from the relative deviation, Fig. 4(a) reveals that ϕ A and ϕ B exchange roles at α ∼ 1.5, around the point where the phase gets its maximum value, see Fig. 3(d).

IV. TOPOLOGICAL PHASE DIAGRAM
The topological state of a system is characterized by the so-called Chern number or topological index [31] with u νk (r) = e −ik·r ψ νk (r) being the periodic part of the Bloch eigenfunctions. Since the band structure of the Haldane model consists on a valence and a conduction band, only the lower energy band enters the sum over occupied states in Eq. (45). In order to efficiently calculate the Chern number, one can rewrite the expression in (45) as where Ω(k) stands for the Berry curvature [32]. This quantity can be accurately computed by means of the Wannier interpolation technique, as discussed in [27,28,33]. In our calculations, we find that a fine 5000 × 5000 k-mesh is required in order to converge the integral of Eq. (46).
The Chern number represents a topological property and takes only integer numbers [31]. Its value is intimately connected to the band structure and the gaps opened by symmetry breaking at the Dirac points. If a gap is opened solely by inversion symmetry breaking, the state of the system is topologically trivial with C = 0. On the other hand, if the gap is opened by time-reversal symmetry breaking, then the system is found in a topologically non-trivial state with C = 0. When both symmetries are broken, the topological state of the system depends on the relative strength of the inversion and time-reversal symmetry breaking.
The topological phase diagram of the Haldane model has been traditionally drawn as a function of ϕ and /|t 1 | [1,18]. In order to facilitate the discussion, let us rewrite here the analytic expression in Eq. (26) that defines the boundary between the different insulating regions, namely In the original formulation, in which the dependence of ϕ on α is derived by means of the Peierls substitution [1,18], the whole phase diagram is accessible. However, since the Peierls substitution is incorrect [20], the possible values of ϕ are actually limited to a finite range that depends on s, as discussed in section III C (see e.g.  Fig. 3(d)). This is shown in Fig. 5, where the accessible region for each value of s is represented by the vertical (dashed) lines. Actually, only a small portion of the nominal phase diagram can be accessed (see the inset), as the maximum allowed values of ϕ are much smaller than π.
In the figure, the dots represent a non-trivial topological state with C = ±1. The fact that almost all these points lie in between the black solid lines [34] proves that -in the allowed accessible region -the phase diagram of the microscopic Hamiltonian is well described by the analytical expression of Eq. (47) for the Haldane model. Owing to the above analysis, we suggest that a more appropriate way to draw the topological phase diagram is in terms of the physical parameters that characterize the underlying continuous Hamiltonian, namely α, χ A and s. This is shown in Fig. 6, where we plot the phase diagram in the α − χ A plane, for three different values of s. Importantly, the Fig. evidences that the topological insulating phase with C = 0 shrinks dramatically as the system becomes more and more tight-binding (that is, by increasing s). Notice that the sign of the Chern number in the topological insulator phase (C = ±1) is consistent with the sign of α, and independent on the sign of χ A . Notice also that the probability of finding the system in the topological insulator phase increases consistently by decreasing the value of |χ A |.
As previously anticipated, the structure of the phase diagram is intimately connected to the behavior of the gaps at the Dirac points. This is illustrated in Fig. 7, where we plot the gaps δ + and δ − as a function of α > 0 and three different values of χ A > 0 for fixed s = 5. Noteworthy, the gap closing does not take place exactly at k − D but in a close-by non-high-symmetry point. The origin of this feature has already been discussed in Sec. II C. For  [35]. We find that the gap closing point is slightly shifted for different values of s, but lies always very close to k − D . In all cases, the deviation from k − D represents a minor correction, and can be safely ignored in the following discussion.
Notably, Fig. 7 reveals that the gap has a maximum at α 1.0 k L , impliying that the effect of the vector potential in opening the gap is limited, as expected from Eq. (26). It is also noteworthy that when χ A is relatively small, as in Figs. 7a and 7b, the gap δ − vanishes for two different values of α (the role of δ + and δ − is exchanged for α < 0) . In fact, owing to the non monotonic behavior of ϕ as a function of α, see Fig. 3d and Ref. [20], there are two different values of α for which Eq. (47) can be satisfied (notice that the two values of ϕ at the phase boundaries may be slightly different due to the fact that t 1 also depends on α). The intermediate region between these two values, which is represented by a grey shaded area in the figures, corresponds to a topological non-trivial state (C = 1) where the effect of time-reversal symmetry breaking is stronger than inversion symmetry breaking. As mentioned above, the smaller the value of χ A , the larger the region with C = 0 as a function of α. By increasing χ A , the topological insulating phase shrinks and eventually disappears, as shown in Fig. 7(c).
To conclude our analysis, let us discuss why the phase diagram of Fig. 6 shrinks as s is increased. For such purpose, in Fig. 8 we illustrate the evolution of the gap as a function of α and s for fixed χ A . This figure evidences that the maximum of the gap decreases as s is increased; in other words, the relative effect of time-reversal symmetry breaking decreases with increasing s. As a consequence, even relatively low values of χ A can avoid gap closing provided s is large enough, as in the case of s = 9 in Fig. 8. This, in turn, implies that the phase transition to the topological insulator phase is restricted to smaller values of χ A as s is increased, in agreement with the phase diagram of Fig. 6.

V. CONCLUSIONS
In summary, we have presented an ab initio analysis of a continuous Hamiltonian [18] that maps into the celebrated Haldane model [1]. The tunnelling coefficients of the tight-binding model have been computed by means of two independent methods, one based on the maximally localized Wannier functions and the other on a closed set of analytical expressions in terms of the energy spectrum at selected high symmetry points in the BZ. The two approaches present a remarkable agreement. In particular, we have shown that the gaps created either by inversion or time-reversal symmetry breaking are very well described by the tight-binding model, which reproduces accurately the exact behavior. In addition, we have calculated the topological phase diagram in terms of the physical parameters entering the microscopic Hamiltonian, finding that only a small portion of the original phase diagram discussed by Haldane can be actually accessed within this model. Moreover, we have shown that the non-trivial topological phase with non-zero Chern number is suppressed as the system enters the deep tight-binding regime. We believe that, besides its conceptual implications, this work is relevant for a possible experimental implementation of the Haldane model following the proposal in Ref. [18]. , as the amplitude α of the vector potential is varied and the system crosses the topological phase boundary. Marzari and Vanderbilt showed that this functional can be divided into three parts, namely Ω = Ω I + Ω D + Ω OD [23]. The term Ω I is gauge-invariant (namely, independent of the choice of the unitary transformations U νν (k) in Eq. (6)), whereas the diagonal term Ω D and the off-diagonal term Ω OD do depend on the gauge choice. In Fig. 9 we show the behavior of the three terms of the spread as a function of α, for fixed values of s and χ A . Here, the non-trivial topological phase is indicated by the grey shaded area. All the components of the spread show a continuous behavior, even across the boundary between the trivial and non-trivial topological states. Then, it is interesting to note that, while the gauge-invariant term Ω I shows a monotonic decrease as a function of α, the gauge-dependent terms Ω D and Ω OD show a non monotonic behavior that is reminiscent of what we observed for the gap (see Fig. 7) and for the complex phase of the next-to-nearest tunnelling coefficient in Fig. 3(d).
We notice that the smooth behavior of the spread shown by our calculations differs from an earlier analysis of MLWFs in the context of the Haldane model performed by Thonhauser and Vanderbilt [36]. There the authors found a breakdown of the usual procedures to build MLWFs as the system approaches the topological phase boundary, resulting in a divergence of the spread functional. The fundamental difference between our approach and the one followed in Ref. 36 resides in the set of bands considered for the construction of the ML-WFs. In fact, whereas our set includes both the valence and conduction bands, their approach included only the valence band. This is a crucial difference, since the net Chern number of a single band in the topological phase is finite, therefore it becomes impossible to choose a smooth periodic k-space gauge of the Bloch orbitals and the procedure for constructing the MLWFs fails. In our case, in contrast, the net sum of the Chern numbers of the valence and conduction bands remains null, hence there is no formal impediment for the construction of the ML-WFs.
requires an energy cutoff of 50 E R , whereas the rest of the terms in the Hamiltonian are converged with 10 E R . Finally, for extracting the tight-binding parameters using the formulas discussed in section III B, we have used a direct diagonalization of H 0 in Eq. (1) by means of a standard Fourier decomposition. In this case, the vector potential term in Eq. (B2) is transformed into a non diagonal matrix in momentum space.