Maximally localized Wannier functions for ultracold atoms in one-dimensional double-well periodic potentials

We discuss a method for constructing generalized Wannier functions that are maximally localized at the minima of a one-dimensional periodic potential with a double-well per unit cell. By following the approach of (Marzari M and Vanderbilt D 1997 Phys. Rev. B 56, 12847), we consider a set of band-mixing Wannier functions with minimal spread, and design a specific two-step gauge transformation of the Bloch functions for a composite two band system. This method is suited to efficiently compute the tight-binding coefficients needed for mapping the continuous system to a discrete lattice model. Their behaviour is analyzed here as a function of the symmetry properties of the double-well (including the possibility of parity-breaking), in a range of feasible experimental parameters.


Introduction
Ultracold atoms in optical lattices are attracting an increasing interest as quantum simulator for condensed matter systems [1,2]. Optical lattices are realized by means of one or more continuous sinusoidal potentials created by exploiting the electric dipole interactions between the atoms and the laser beams [1]. Depending on the beam geometry one can realize one-, two-, or three-dimensional periodic lattices, with one or more wells per unit cell [3]; quasiperiodic structures are also possible (see e.g. [2,4,5,6,7]).
In the tight binding regime, when the lattice intensity is sufficiently high so that the atoms are deeply localized in the lowest vibrational states of the potential wells (each well being associated to a site of a discrete lattice), it is convenient to map the system hamiltonian onto a discrete lattice model (or tight binding model ), representing the usual tool for theoretical calculations. The paradigms are the Hubbard model for fermions [8], and the Bose-Hubbard model for bosons [9]. These models are characterized by tunnelling coefficients related to the hopping between neighbouring sites, and interaction strengths which characterize the onsite interaction among the atoms, whose actual values depend on the parameters of the underlying continuous model.
As the mapping between the continuous and discrete versions of the system hamiltonian is achieved by means of an expansion over a basis of localized functions at each potential well [1,10,11], a precise knowledge of these basis functions is therefore important to connect the actual experimental parameters with the coefficients of the discrete model employed in the theoretical calculations.
For a simple sinusoidal potential, the natural basis is provided by the exponentially decaying Wannier functions discussed by Kohn [12,13]. Notably, in this case the expression for the tunnelling coefficient turns out to depend just on the Bloch spectrum [14], being therefore independent on the basis choice (this does not hold for the interaction coupling). Analytic expressions for both coefficients can be obtained by means of different approximations [11,15].
In the case of two wells per unit cell, the Kohn-Wannier recipe is not sufficient. For example, for a symmetric double-well the Kohn-Wannier functions display the same symmetry as the local potential structure [13,16] and thus they occupy both wells and cannot be associated to a single lattice site. A common approach used in the literature is that of the so-called atomic orbitals [19,17,18], that has been recently employed e.g. for the case of a symmetric double-well unit cell of two-dimensional graphene-like optical lattices [20]. This method is based on a specific ansatz, according to which tight-binding Wannier functions are constructed from linear combinations of wave functions deeply localized in the two potential wells of the unit cell.
A more general approach is the one proposed by Marzari and Vanderbilt [21], where maximally localized Wannier functions (MLWFs) are obtained by minimizing the spread of a set of Wannier functions by means of a suitable gauge transformation of the Bloch eigenfunctions. This method coincides with the Kohn method for the single band MLWFs of a one-dimensional potential, but it can be extended to more complex situations when generalized MLWFs for composite bands are needed. The method is implemented by means of a software package, and is largely employed for computing MLWFs of real condensed matter systems [22].
In this article, we consider the case of a one-dimensional periodic potential with a double-well per unit cell [23,24,25,26,27,28,29], discussing an alternative method for constructing the low-lying generalized MLWFs based on the minimal spread requirement of Marzari et al. [21], specifically suited for double-well potentials. In particular, we consider a composite band formed by the two lowest Bloch bands and, differently from [21], we design a two-step gauge transformation specific for a composite two band system, that can be solved by integrating a set of ordinary differential equations, with suitable boundary conditions. This allows to efficiently compute the tunnelling coefficients and the other tight binding coefficients in terms of the parameter of the continuous potential. We also remark that, though the approach of Marzari et al. [21] was proposed for degenerate bands, we find that the method works properly in a wider regime, provided that the first band gap does not exceed the distance between the second and third band.
The paper is organized as follows. In Section 2 we review the mapping of a manybody hamiltonian onto a tight binding model, discussing in particular the case of a sinusoidal periodic potential with a double well in the unit cell. Then, in Section 3 we describe the gauge transformation for constructing the generalized MLWFs, giving some examples and comparing them with those of the single band approach. Then, in Section 4, we discuss the regime of validity of the composite band approach, and compare the predictions of the full and nearest-neighbour versions of the model by discussing the behaviour of the tunnelling coefficients and the other tight-binding parameters. Final considerations are drawn in the Conclusions. Technical points regarding the mapping on momentum space and the numerical implementation are discussed in the Appendices.

Tight binding models for double well optical lattices
Let us start by reviewing how tight-binding models are defined from a continuous potential. As a specific example, here we consider a one-dimensional many-body hamiltonian for ultracold bosons [10] withψ(x) being the bosonic field operator,Ĥ 0 = −( 2 /2m)∇ 2 + V (x) the single particle hamiltonian, and V (x) a double-well periodic potential of the form with k B = π/d and V 1 strictly non vanishing (V 1 > 0) in order to fix the overall period to d, V (x+d) = V (x). This is a typical potential used in the experiments with ultracold atoms [1], by which one can construct an optical lattice with one or two wells in the unit cell. In particular, the condition for having two minima in each period for any value of θ 0 is V 2 > 0.5V 1 . In the following, the potential amplitudes V i will be expressed in units of E R = 2 k 2 B /2m, the so-called recoil energy for an atom absorbing a photon of the first lattice.
The phases θ 0 and φ 0 are arbitrary; φ 0 represents a rigid shift of the whole potential, while θ 0 can be varied along with the ratio of the amplitudes V 2 /V 1 to change the landscape of the potential. For convenience, the unit cell is defined as having two maxima at the cell borders, with both minima inside the cell (see figure 1). In addition, the angle φ 0 is tuned in order to have a unit cell centered in x = 0, x ∈ [−d/2, d/2], with the absolute minimum in the left well ‡. Depending on the value of θ 0 , it is possible to realize three different configurations, as shown in figure 1. (a) A unit cell with two different minima and with degenerate maxima, for θ 0 = nπ (n ∈ Z); in this case the whole periodic potential has two (classes of) parity centers, corresponding to the two Figure 1. The three possible configurations for the unit cell of the potential in (2) (here V 2 = 2V 1 ): (a) two different minima, with the overall potential having two centers of parity, for θ 0 = 0 (φ 0 ≃ π/4); (b) an asymmetric double-well with parity that is broken globally -in this example θ 0 = π/4 (φ 0 ≃ π/8); (c) a symmetric double-well, θ 0 = π/2 (φ 0 = 0). The black dots in (a), (c) represent the parity centers of the whole periodic potential.
minima, with all the maxima being degenerate. (b) An asymmetric double-well with parity that is globally broken, for any θ 0 ∈ (0, π/2) + nπ/2. (c) A symmetric double-well in the unit cell, for θ 0 = π/2 + nπ; in this case the potential has again two centers of parity, now at the two maxima, with all the minima being degenerate. As anticipated in the Introduction, when the potential wells are deep enough, it may be convenient to map the hamiltonian (1) onto a tight binding model on the discrete lattice corresponding to the potential minima, by expanding the field operatorψ(x) on a basis {f nj (x)} of functions localized around each minimum whereâ † nj (â nj ) represent the creation (destruction) operator of a single particle at site j, and satisfy the usual commutation rules [â nj ,â † n ′ j ′ ] = δ jj ′ δ nn ′ (following from those for the fieldψ).
In the presence of a single well per unit cell, it is known that a basis of localized functions is provided by the exponentially decaying Wannier functions w nj (x) discussed by Kohn [13,14]. In general, this is not the case when there are two wells per unit cell. For example, for a symmetric double-well the Kohn-Wannier functions display the same symmetry of the local potential structure [13,16] and thus they cannot be associated to a single lattice site as as they occupy both wells in the unit cell. In the next section we will show that when the two lowest Bloch bands are sufficiently close to each other with respect to the third band (as it will be clear from the discussion in Sect. 3.2 and 4), we can construct a set of generalized Wannier functionsw nj (x) that are maximally localized at each minima, by following the approach of Marzari and Vanderbilt [21] for a composite band. This corresponds to the generalization of the single band approximation (in case of a single well lattice) to the double well case, as we need at least two localized functions in each lattice cell to map the system on the discrete lattice. Then, in Section 4, we will explore the range of validity of this composite band approaches, highlighting the different implications on the structure of different tight binding models. In the following we will restrict the analysis to the two lowest energy bands, in analogy with the single band approximation for the Bose-Hubbard model [11]. Then, within this approximation, the single particle hamiltonian can be written aŝ where j is the unit cell index whereas ν = A, B substitutes the band index n = 1, 2 being an internal index labelling the left and right sub-wells respectively (see figure 2).
Here the expansion coefficients correspond to the onsite energies E ν = f jν |Ĥ 0 |f jν , and to the tunnelling amplitudes between different (sub)wells T jj ′ νν ′ ≡ − f jν |Ĥ 0 |f j ′ ν ′ . In general, it is customary to further approximate the above expression by neglecting the coupling beyond nearest neighbours both for the single-well [10] and double-well lattices [25,28,29]. This is a reasonable assumption for a single well lattice in the tight binding regime [14,34], but may not be fully justified in the range of the typical experimental parameters for a double well, as we will discuss later on. For this reason, here we keep all the terms corresponding to nearest-neighbouring cells, characterized by the following tunnelling coefficients as shown in figure 2. The corresponding hamiltonian iŝ When the tunnelling J ν and J AB + are negligible the above expression can be further simplified by retaining just the coupling between nearest neighbouring wellŝ This is the analogous of the nearest-neighbours approximation for the single-well case and it is commonly used in the literature [28,29]. Hereinafter, we will refer to the approximate hamiltonians in (8) and (9) as the (single particle) full and nearestneighbour tight-binding models, respectively. Finally, the general form of the interaction term iŝ Here we consider just the onsite contributions in the two sub-wells A and B with U ν = g dx |f jν (x)| 4 . This corresponds to the usual approximation used for the single well case (see e.g. [10]), and is justified when the basis functions {f nj (x)} are strongly localized in each sub-well. In general, once these functions are known, the next-to-leading interaction terms can be straightforwardly computed by evaluating the corresponding superposition integral of four f nj 's.

Single particle spectrum
Let us now consider the single particle spectrum of the hamiltonian (8). By defininĝ the full tight binding hamiltonian (8) can be rewritten aŝ and where the operatorsb kν satisfy the canonical commutation relations, Then, by diagonalizing the matrix h(k), and defining ǫ ± (k) ≡ (ǫ A (k) ± ǫ B (k))/2, we get (see also [20]) that represents the spectrum of the full tight binding model in (8). In addition, with J ν = 0 = J AB+ , the same expression gives also the spectrum for the nearest-neighbour approximation in (9). Instead, in the single band case we simply have [1] ε sb n (k) = E sb n − 2J sb n cos(kd) with ε n (k) being the exact Bloch spectrum; notably these expressions do not depend on the choice of the Wannier basis.

Generalized Wannier functions
In this section we discuss the method for constructing the MLWFs for the double well case. In order to fix the notations, let us first recall some basic properties of periodic systems [19,30]. Owing to the Bloch's theorem, the eigenfunctions of the single particle hamiltonianĤ 0 can be written as ψ nk (x) = e ikx u nk (x), the u nk (x)'s having the same periodicity of the potential and satisfying the following normalization in the unit cell, u mk |u nk = (d/2π)δ mn . The Wannier functions for a single band are defined as with B indicating the first Brillouin zone, k ∈ [−k B , k B ], and R j ≡ jd, whereas generalized Wannier functions for composite bands have the same formal expression but are built up from a linear combination of Bloch eigenstates, namelỹ with UU † = 1 and U nm (k + 2k B ) = U nm (k). The Wannier functions satisfy the ortho-normality relation w nj |w n ′ j ′ = w nj |w n ′ j ′ = δ nn ′ δ jj ′ . We also remark that the generalized Bloch functionsψ nk are not eigenstates ofĤ 0 ; however, their ortho-normality relations are preserved, owing to the unitarity of the transformation matrices U nm (k). Different choices of the matrices U nm (k) lead to different results and it is customary to speak about gauge dependence of the Wannier functions due to the k-dependence of the U(n) transformations. In the single band case the arbitrariness reduces to the abelian U(1) group of phase transformations coming from the freedom in choosing the Bloch basis at each k.
A general approach to obtain MLWFs has been proposed in a seminal paper by Marzari and Vanderbilt [21], where MLWFs are obtained by minimizing the generalized Wannier spread Ω = n [ x 2 n − x 2 n ] by means of a suitable gauge transformation of the Bloch eigenfunctions. In the case of a single band the method returns the Kohn result [21]. Marzari and Vanderbilt have shown that the spread can be written as the sum of two positive terms, Ω = Ω I +Ω, the first being gauge invariant and therefore fixing the minimal spread. The gauge dependent termΩ can be further split into the diagonal and off-diagonal components,Ω = Ω D + Ω OD , that in the one-dimensional case read Both Ω D and Ω OD can be written in terms of the generalized Berry vector potentials A nm (k), defined as [32,33] We also recall that the integral over the first Brillouin zone of A nn (k) gives the one band Zak-Berry phases γ n which are proportional to the offset of the Wannier function centers, [31,21], yielding x n0 = A nn B (this relation is preserved under a generic unitary gauge transformation, as in (22)). It is also worth to remember that the single Wannier centers are invariant only under single band U(1) gauge transformations, whereas for general U(n) transformations only their sum is conserved [21]. Then, the expressions for Ω D and Ω OD can be written as and in general can be reduced by means of a functional minimization in k space, as discussed in [21,22].
Here we use a different approach, specifically suited for constructing the set of MLWFs for the double-well case. Namely, we show that the minimization problem can be reformulated by identifying a specific gauge transformation for a composite band in one dimension, expressed in terms of a set of ordinary differential equations with periodic boundary conditions. We recall that in one dimensionΩ can be made strictly vanishing, and this corresponds to find a gauge (also called parallel transport gauge) in which the matrix A nm (k) is diagonal, with the diagonal elements being constant and equal to their mean values. The latter are related to the eigenvalues of the matrix generalizing the Berry phase to the non abelian case [21].

Gauge transformations and differential equations
The diagonal and off-diagonal spreads Ω D and Ω OD can be minimized either simultaneously or independently. For the following discussion, it is useful to distinguish between two kind of gauge transformations.
I. Diagonal U(n) transformations which correspond to a set of single band gauge transformations of the form |u nk → |ũ nk = e iφn(k) |u nk (29) with φ n (k) being a real continuous (differentiable) function of k, such that φ n (k + 2k B ) = φ n (k) + 2πℓ (ℓ integer) in order to have periodic and single valued Bloch eigenstates. Then, we can set ℓ = 0 without loss of generality. As discussed in [21], this transformation may be used to minimize each term of Ω D , as it affects the Wannier spread Ω n = x 2 n − x 2 n while preserving their centers x n (modulo a lattice vector). In particular, Ω D can be set exactly to zero [21]. In fact, we have A nn (k) → A nn (k) − ∂ k φ n (k), and therefore that vanishes by imposing This equation can be readily solved numerically, as discussed in Appendix A and Appendix B. In addition, it is straightforward to verify that under the transformation (29), A 12 (k) changes only by a phase factor, and therefore Ω OD remains unchanged. II. A full gauge transformation in the composite band of the form where, as already said, U(k) ∈ U(N) (for a N-composite band), constrained to the following periodic condition Under such a transformation the generalized Berry potentials transform as In general, U(N) can be written as a semidirect product SU(N) ⋊ U(1), with U(1) being subgroup of U(N) consisting of matrices of the form diag(1, 1, . . . , e iχ ). As anticipated, here we restrict to a 2 × 2 case, and therefore we have with |z 1 | 2 + |z 3 | 2 = 1, r(k) = e iχ(k) . Moreover, by using the following parametrization for a matrix S ∈ SU(2), S = e iα σ ·n/2 withn = (cos ϕ sin θ, sin ϕ sin θ, cos θ) and σ i being the Pauli matrices, we can write with χ = χ(k), ϕ = ϕ(k), α = α(k) and θ = θ(k).
Then, since Ω transforms as follows in order to getΩ = 0 one has to impose (see (34)) Notice that the former, (39), is an integro-differential equation where the right-end side corresponds to the center of the Wannier functions, Ã nn B = w nj |x|w nj , that are not known a priori (only their sum is conserved in the parallel transport gauge). Therefore, in the following we will consider a specific transformation that makes Ω OD vanish without specific requirements on the transformation of Ω D , given by the solution of (40). Then, by using (36), the latter can be transformed in a system of four differential equations for α, θ, ϕ, χ, whose normal form is ∂ k θ = cos θ sin α sin 2 (α/2) (A R 12 cos η + A I 12 sin η) where we have defined η ≡ ϕ − χ, with ∂ k η = 0. The solution of (43) is χ = χ 0 , ϕ = ϕ 0 . Then, it is evident that only two equations are left, namely (41) and (42), with η = ϕ 0 − χ 0 playing the role of a parameter, and with one of the angles among ϕ 0 and χ 0 being arbitrary. We then can choose χ 0 = 0 without loss of generality. In addition, in order to conform to (33), the angles α and θ must satisfy the following periodicity conditions with ℓ, ℓ ′ integers which can be taken ℓ = ℓ ′ = 0 without loss of generality. This set of equations, (41), (42), (44), and (45) can be solved as discussed in Appendix B. The resulting transformation is an SU(2) matrix S(k) of the form (36) with χ = 0 and ϕ constant. Summarizing, the procedure to makeΩ vanish can be divided in two steps: (a) a gauge transformation of type II to make Ω OD vanish (without specific requirements on the transformation of the diagonal elements A nn (k)); (b) a set of two gauge transformation of type I (one for each band) that makes Ω D vanish without affecting Ω OD . Therefore, the full transformation can be decomposed in the following way It is straightforward to verify that (46) can be cast again in the form (35) by properly redefining the parameters z 1 , z 3 , r.

Examples of MLWFs
Here we present some example of generalized MLWFs obtained with the transformation discussed above, and we compare them with the single band MLWFs. The latter can be obtained by using just the gauge transformation of type I, and correspond to the exponentially decaying Wannier functions discussed by Kohn [13,21]. Both gauge transformations are solved by using the representation of Bloch functions in k-space, by means of the numerical methods discussed in Appendix A and Appendix B.
In figure 3 we show the case of a symmetric double well, θ 0 = π/2. In this case the single band MLWFs have the same symmetry of the potential [13], therefore occupying both wells in the unit cell. Instead, each of the generalized MLWFs nicely localizes in one of the sub-wells, as wanted §. § We remark that in this symmetric case one could construct a set of localized functions resembling the generalized MLWFs by considering symmetric and antisymmetric combinations of the single band MLWFs. These combinations (not shown) nicely reproduce the bulk density of the generalized MLWFs when each unit cell can be regarded as a single double-well (that is, large barriers at the cell borders). However, this approximation does not work in the region of the tails (those of the generalized MLWFs decaying much faster), therefore not being reliable for computing the tunnelling coefficients. It is also   Then, in figures 4, 5 we show two cases of an asymmetric double well, here for θ 0 = 0. These figures are useful for illustrating the role of the band gaps. The case in figure 4 corresponds to a band gap between the first two bands that is larger than the one between the second and third band; in this case both the single band MLWFs are already localized within a single well, the one in panels (b,d) having a small residual amplitude around the neighbouring wells. The generalized MLWFs do not differ much from the former in linear scale, the effect of the gauge mixing transformation being a reduction of the lateral lobes of the single band MLWF in (b,d), but the price to pay is a substantial increase of the width of the other one in log scale, see panel (a). Actually, in this case composite band approach is not fully justified, due to the large band gap possible to prove analytically that this approximation implies J AB− = J AB+ , that is manifestly in contradiction with the definition of the tunneling coefficients given in figure 2. Nevertheless, such wave functions may be useful for discussing the boundary conditions for (41)-(43), see Appendix B. between the first and second band, and one should consider the structure of the upper bands. This however goes beyond the scope of this work. Instead, figure 5 shows a case still for θ 0 = 0 but where an almost degeneracy between the two minima (that corresponds to a quasi degeneracy of the lowest two bands) is restored thanks to a lower value of V 1 (here V 1 = 1 instead of 10 as in the former figure). Here the advantage of using the generalized MLWFs is clearly visible.

Tight binding regimes and tunnelling coefficients
In this section we will discuss the features of the composite band approach, by comparing the predictions of the full and nearest-neighbour versions of the model, and discussing the behaviour of the tunnelling coefficients and the other tight-binding parameters.
We recall that the potential is characterized by three parameters, V 1 /E R , V 2 /E R and θ 0 . The latter can be restricted to the range [0, π/2] without loss of generality, see figure 1. As for the potential amplitudes, we chose 10 ≤ V 2 /E R ≤ 40 in order to fulfil the tight binding regime, in a feasible range for current experiments with ultra cold atoms [1]. In addition, for having both wells deep enough in order to have the corresponding MLWFs mostly localized within each well, V 1 should not exceed V 2 /2.
As we have anticipated in the previous section, in principle the composite band approach is justified in a situation of quasi degeneracy between the first two Bloch bands, that is when their separation ε g 12 is much smaller that the band-gap ε g 23 between the second and third band. Therefore, it is convenient to define the ratio R ≡ ε g 12 /ε g 23 , whose behaviour is shown as a density plot in figure 6 as a function of θ 0 and V 2 , for different values of V 1 . As a matter of fact, we find that the composite band approach provides a basis of functions well localized in each of the two sub-wells and correctly reproduces the lowest two Bloch energy bands even up to R ≈ 1, in the parameter range discussed above. In order to characterize the level of fidelity in reproducing the exact single particle Bloch spectrum (that can be readily computed as discussed in Appendix A), we define the following quantity that represents the ratio of the quadratic spread between the exact Bloch spectrum ε n (k) and that of the tight binding hamiltonians (8) and (9), to the bandwidth ∆ε n ≡ (ε max n − ε min n ). We also notice that the formula (17) for ε tb n (k) provides the same numerical result of the expression (C.6) for the energy bands in terms of the gauge transformations, and that these formulas have a better accuracy in reproducing the exact Bloch spectrum than the single band expression (17), especially in the region close to the symmetric case θ 0 = π/2. The quantity δε n is shown in figure 7 for V 1 = 5, as a function of θ 0 (V 2 = 20) and V 2 (θ 0 = π/2), in the left and right panels respectively. These figures refer to a horizontal and a vertical cut in the left panel of figure 6, and confirm that in the regime R 1 the full tight binding model reproduces the exact energies with great accuracy. As for the nearest neighbour approximation commonly used in the literature, in general this is not the case, as it works reasonably only for R 0.1, that is for θ 0 ≃ π/2 and large V 2 . Therefore, attention must be payed to the regimes where it is allowed to neglect some of the next-to-nearest tunnelling terms.

Tunnelling coefficients
Let us now consider the explicit expressions for the onsite energies and the tunnelling coefficients in (8). They can be expressed in terms of the gauge transformations discussed in the previous section as (see (21) and (46))  Figure 9. Plot of the modulus of T AB (×100) and of the onsite energy difference E B − E A , for the same parameters of the previous figure.
with ν = A, B = 1, 2 and ∆φ(k) = φ 2 (k) − φ 1 (k). We recall that all these terms are relevant for the full tight-binding model, while the nearest-neighbour approximation commonly used in the literature corresponds to retain just the contribution of E ν , T AB and J AB − . We also notice that all the above parameters are gauge dependent; E ν and J ν depend just on the gauge mixing transformation; T AB and J AB ± on the whole gauge transformation.
It is worth mentioning that in the parallel transport gauge only the modulus of the tunnelling coefficients is univocally defined. However, as it is briefly discussed in Appendix B, there is the freedom to choose them real.
The behaviour of the tunnelling coefficients (in modulus, rescaled to T AB ) is shown in figure 8 as a function of θ 0 and V 2 (for V 2 = 20 and θ 0 = π/2, respectively; here  Figure 10. Plot of the rescaled interaction term I ν , for V 1 = 5, as a function of θ 0 for V 2 = 20 (a), and as a function of V 2 for θ 0 = π/2 (b). Empty squares: first band; Solid squares: second band. In the right panel the two lines are superimposed owing to the double well symmetry.
V 1 = 5), whereas the absolute variation of T AB and of the onsite energy difference E B − E A is shown in figure 9. The former figure reveals that the dependence on θ 0 is rather weak, the only notable effect being that at θ 0 = 0 (where all the maxima are degenerate) T AB = J AB+ , whereas at θ 0 = π/2 we have J A = J B (and E A = E B , as in this case the minima are degenerate). Instead, the tunnelling coefficients J A , J B and J AB+ are substantially suppressed with respect to T AB and J AB− as V 2 is increased, and consequently the results of the nearest-neighbour approximation improve. Also, figure  7a shows that the bands energies obtained with the nearest-neighbour approximation better reproduce the exact results by approaching θ 0 = π/2 in spite of the very weak variation of all the tunnelling coefficients shown in figure 8. Actually, a detailed analysis reveals that these results are due to the reduction of the band gap ε g 12 and to the increase of the lowest bands widths ∆ε 1,2 when going from θ 0 = 0 towards θ 0 = π/2. Both from (17) and (C.6) it is possible to realize that the relative weight of the terms containing J A and J B becomes consequently less relevant. Finally, in figure 10 we consider the behaviour of the rescaled interaction term (see (11)), I ν ≡ U ν /g = dx |w jν (x)| 4 , again as a function of θ 0 and V 2 , for the parameters of the preceding figures. This term is also called inverse participation ratio (IPR), and represent the number of occupied lattice sites in the case of a model defined on a discrete lattice (see e.g. [7,36]). In the present continuous case, values of I ν close or greater than 1 imply that the MLWFs are localized on a distance smaller than the lattice spacing .

Conclusions
We have presented a method for constructing a set of maximally localized Wannier functions (MLWFs) for a one-dimensional periodic potential with a double-well per unit cell. Starting from the minimal spread requirement by Marzari et al. [21], In the regime of parameter considered here, the next-to-leading interaction terms -see (10) -are suppressed by about 2 orders of magnitude.
we have designed a two-step gauge transformation specific for a composite two band system, consisting a set of ordinary differential equations with periodic boundary conditions. In the tight binding regime, this method provides a proper set of MLWFs when both wells are sufficiently deep, provided that the band gap between the first two bands is larger than the one between the second and third band. By using these MLWFs one can map the continuous hamiltonian onto a tight binding model, whose coefficients are expressed in terms of the above gauge transformations, providing a direct way to relate the tight binding coefficients to the experimental parameters characterizing the continuous potential. We have also shown that in general, in the range of typical experimental parameters, it is necessary to retain all possible tunnelling couplings between the different sub-wells of neighbouring cells (full model), whereas the usual nearest neighbouring approximation used in the literature is justified only for V 2 /E R , V 2 /V 1 ≫ 1 and θ 0 ≃ π/2. The method can be applied also in higher dimensions, in case of separable lattice potentials.
Since the diagonalization at each point in k-space is uncorrelated from that at neighbouring ks, in general the diagonalization routines return a set of coefficientsc mℓ (k) with an arbitrary phase factor e iϕ nk that cannot be controlled a priori, and that affects their differentiability in (A.3). However, a set of smooth (differentiable) c mℓ (k) can be obtained as follows. We notice that in x = 0 the periodicity condition (A.4) reads u nk+2k B (0) = u k (0). Therefore, given a set of Bloch eigenfunctionsū nk (x) generated numerically (affected by the same spurious phase factors e iϕ nk of thec nℓ (k), see (A.1)) it is possible to define a new set of smooth Bloch eigenfunctions u nk (x) by fixing all the phases to zero at x = 0 (i.e. by imposing them to be real in x = 0), therefore fixing the spurious phases. Formally, this corresponds to a diagonal gauge transformation generated by φ n (k) = − arg(ū nk (0)) = − arg ( ℓc nℓ (k)). Since φ n (k) does not depend on ℓ, the c nℓ (k) transform in the same way and this can be used as a prescription for ensuring the smoothness of the c nℓ (k) and u nk as a function of k.
Gauge mixing transformation: method 2. In order to solve (41), (42), we have also employed another method that originates from the observation that in the case of two degenerate symmetric minima in the elementary cell (θ 0 = π/2), simple symmetric and antisymmetric combinations of Bloch eigenstates give rise to two Wannier functions located each at one of the two degenerate minima (though they are not the maximally localized ones). In this situation, we have verified that initial conditions for the matrix S nm (±k B ) -corresponding to symmetric and antisymmetric combinations of Bloch eigenstates -allow to obtain the periodic solutions minimizing Ω OD which coincide both integrating from −k B to k B and vice-versa. By smoothly changing θ 0 , we start from the same initial conditions and then we iterate the integration of (41), (42) by taking as initial angles, at each step, the mean values of the final values of the preceding integration (at the two boundaries of the BZ). The procedure rapidly converges with Ω OD reduced by several order of magnitude (this depends on the number of grid points and on the region of parameters, as for the previous method). In particular, for θ 0 = π/2 (symmetric double well) we have always found solutions with ϕ 0 = π/2 and the following boundary conditions for α and θ α(±k B ) = θ(±k B ) = π/2. (B.1) By smoothly changing θ 0 , the corresponding boundary conditions are smoothly connected with this set (in the range of V 1 and V 2 that we have explored). Notably, this choice yields real tunnelling coefficients (in principle this is not guaranteed, see Appendix C.2).

Appendix C. Uniqueness of the physical results
Appendix C.1. Energy bands in the truncated expansion Let us start by considering the following exact expression for the Bloch spectrum ε n (k) = ℓ e ikℓd T n (ℓ) (C.1) with the coefficients T n (ℓ) being the expectation values ofĤ 0 on single band Wannier functions [14] T n (ℓ) ≡ w nj |Ĥ 0 |w nj+ℓ = d 2π B dk ε n (k)e −ikdℓ .

(C.2)
This expression is manifestly gauge invariant. Instead, the analogous expression in terms of the generalized MLWFs depends on the gauge. In fact, with some manipulations we can write (here we are implicitly considering just the first two bands, n = 1, 2) ε n (k) = ℓ e ikℓdP n (k, ℓ) (C.3) withP n (k, ℓ) = n ′ n ′′ U n ′′ n (k)U * n ′ n (k)T n ′′ n ′ (ℓ) (C.4) andT n ′′ n ′ (ℓ) = w n ′′ j |Ĥ 0 |w n ′ j+ℓ . (C.5) By truncating (C.3) to nearest neighbouring cells (ℓ = 0, ±1) and using the definitions of the onsite energies and of the tunnelling coefficients in (5)- (7), we have with U nm (k) = e iφn(k) S nm (k) and S ∈ SU(2) (see Sect. 3). Though each energy ε n (k) in (C.6) depends on the gauge (through the the matrix U nm (k)), it is not difficult to demonstrate that their sum, n=1,2 ε n (k), is gauge invariant like the analogous expression obtained form (C.1). In addition, once we have chosen the parallel transport gauge, corresponding to the generalized MLWFs with minimal spread (Ω = 0), the expression (C.6) is univocally defined despite the presence of residual phase arbitrariness (see next paragraph).

Appendix C.2. Residual phase ambiguities
Let us now show that the relevant quantities are univocally defined in the parallel transport gauge, despite possible ambiguities in the procedure for obtaining the matrix U nm (k). First of all, we remind that the whole procedure starts from a particular gauge, as discussed in Appendix A. Choosing a different initial gauge amounts to perform a diagonal transformation on Bloch eigenstates |ψ mk → e iφm(k) |ψ mk (C.7) being φ m (k) an arbitrary periodic real function (see Sect. 3). Consequently, the full gauge transformation matrix leading to the parallel transport gauge changes as U nm (k) → U nm (k)e −iφm(k) . (C.8) Then, it is easy to see that the onsite energies and tunnelling coefficients in (5)-(7) are invariant under this transformation, as well as the energies in (C.6), whereas the real and imaginary parts of the MLWFs change, conserving their modulus. The procedure to makeΩ vanishing does not imply any other k-dependent ambiguity on the matrix U nm (k). In addition, in the parallel transport gauge we still have the freedom to perform two independent U(1) gauge transformations for each band with angles λ n which necessarily obey ∂λ n /∂k = 0. These terms appear as multiplicative factors on the rows of the matrix U nm . Finally, it turns out that -given a set of parameters V 1 , V 2 and θ 0 -there is not a unique set of boundary conditions for the angles ϕ 0 , α and θ satisfying (41), (42) (and thus not an unique transformation S nm (±k B ) leading to Ω OD = 0). For θ 0 = π/2 these boundary values can be fixed as discussed in Appendix B. Different choices may also be possible. However, it is not difficult to prove that none of these ambiguities affect either the energies in (C.6), or the modulus of the MLWFs and of the tunnelling coefficients. Instead, their real and imaginary parts are not univocally defined and can be chosen somewhat arbitrarily. In particular, real tunnelling coefficients can be obtained by using the set of boundary conditions given at the end of Appendix B.