Coherent superposition of current flows in an atomtronic quantum interference device

The progress achieved in optical micro-fabrication has led to the foundation of atomtronics: Bose–Einstein condensates manipulated with lithographic precision in optical circuits of various intensities and spatial shapes [1–5]. The neutrality of the atoms carrying the current (substantially reducing decoherence sources), the flexibility on their statistics (bosonic/fermionic) and interactions (tunable from short to long-range, from attractive to repulsive) are some of the key features of atomtronic circuits. Atomtronics sets a new stage for quantum simulations [6], with remarkable spin-offs in other fields of science and technology. This activity is believed to lead, in turn, to an improved understanding of actual electronic systems.

An important representative example of an atomtronic circuit is provided by a Bose–Einstein condensate flowing in a ring-shaped trapping potential [7–15]. A barrier potential painted along the ring originates a weak link, acting as a source of back-scattering for the propagating condensate, thus creating an interference state with the forward scattered current. This gives rise to an atomic condensate counterpart of the celebrated rf-SQUID—a superconducting ring interrupted by a Josephson junction [16, 17], namely an atomtronic quantum interference device (AQUID). Due to the promising combination of advantages characterizing Josephson junctions and cold atoms, the AQUID has been the object of recent investigations [18, 19]. The first experimental realizations made use of a Bose–Einstein condensate free to move along a toroidal potential, except through a small spatial region, where an effective potential constriction (giving rise to the aforementioned weak link) is created via a very focused blue-detuned laser or via a painting potential [20–23].

On the theoretical side, it has been demonstrated that the two currents flowing in the AQUID can, indeed, define an effective two-level system, that is, the cold-atom analog of flux qubits [24, 25]. The potential constriction breaks the Galilean invariance and splits the qubit levels that otherwise would be perfectly degenerate at half-flux quantum. In this context, it is of vital importance for the qubit dynamics that a good energy resolution of the two levels could be achieved in realistic physical situations (while keeping the qubit well separated from the rest of the many-body spectrum).

In this paper we focus on ring-shaped confinements with a lattice modulation and a potential constriction. This set-up, that can be realized following different routes (see, e.g., [25]), presents several advantages for the design of an AQUID. First of all, assuming that the bosons occupy only the lowest Bloch band ⁸ , the ring lattice helps in controlling the current. For instance, because of the one-dimensional dynamics, the vortex formation rate along the flow is negligible. Secondly, it helps localizing the barrier effect to a point-like scale with respect to lattice spacing, which should in turn yield a favorable scaling of the qubit gap with the bosonic density [27]. Moreover, it provides an easy route to realize interacting ring–ring architectures [14, 25] ⁹ .

This issue has been considered so far only in some limiting cases, e.g. for particular types of superposition states or in the infinitely strong interacting regime [27, 29]. We perform a systematic study on the quality of the qubit in the cold-atom ring lattice: in particular, we characterize the energy structure at the degeneracy point at half-flux quantum, and study how it is possible to observe experimentally the superposition of the current flows. By employing a combination of analytical and numerical techniques, that allows us to cover all the relevant physical regimes of system sizes, filling, barrier and interaction strengths, we show that: (i) the gap $\Delta {{E}_{1}}$ between the states of the effective two-level system scales as a power law with the system size; (ii) at a mesoscopic scale, a qubit is well-defined, with $\Delta {{E}_{1}}$ displaying a favorable dependence in a wide range of system parameters; (iii) the superposition state is detectable in the momentum distribution of the bosonic gas, which is measurable via time-of-flight (TOF) expansion, and (iv) the momentum distribution exhibits a subtle interplay between barrier strength and interaction.

The paper is organized as follows. In the next section we present the physical system of interacting bosons on a 1D ring lattice with a potential constriction, and the effective two-level system giving rise to the AQUID. In section 3 the energy spectrum of the system and its scaling with system size, filling, and interaction is analyzed. In section 4 we show how the state of the AQUID can be read out through TOF expansion images of the gas. Finally, we draw our conclusions in section 5. Technicalities on the employed methods and further details are provided in the appendices.

We consider a system of N interacting bosons at zero temperature, loaded into a 1D ring-shaped optical lattice of M sites. The discrete rotational symmetry of the lattice ring is broken by the presence of a localized potential on one lattice site, which gives rise to a weak link. The ring is pierced by an artificial (dimensionless) magnetic flux $\Omega$, which can be experimentally induced for neutral atoms as a Coriolis flux by rotating the lattice at constant velocity [22, 30], or as a synthetic gauge flux by imparting a geometric phase directly to the atoms via suitably designed laser fields [31–33].

In the tight-binding approximation, this system is described by the 1D Bose–Hubbard (BH) Hamiltonian

$$\begin{eqnarray}&&H=\mathop{\sum }\limits_{j=1}^{M}\left[ -t({{{\rm e}}^{-{\rm i}\Omega /M}}b_{j}^{\dagger }{{b}_{j+1}}+{\rm h}.{\rm c}.)+\frac{U}{2}{{n}_{j}}({{n}_{j}}-1)+{{\Lambda }_{j}}{{n}_{j}} \right],\end{eqnarray} \tag{ 1 }$$

where ${{b}_{j}}\;(b_{j}^{\dagger })$ are bosonic annihilation (creation) operators on the jth site and ${{n}_{j}}=b_{j}^{\dagger }{{b}_{j}}$ is the corresponding number operator. Periodic boundaries are imposed by taking ${{b}_{M+1}}\equiv {{b}_{1}}$. The parameter U takes into account the finite scattering length for the atomic two-body collisions on the same site; ${{\Lambda }_{j}}$ defines an externally applied local potential. Periodic boundary conditions are assumed in order to account for the multiply connected geometry of the ring system. The presence of the flux $\Omega$ is taken into account through the Peierls substitution: $t\to t{{{\rm e}}^{-{\rm i}\Omega /M}}$ (t is the hopping amplitude). In the thermodynamic limit, the BH model for ${{\Lambda }_{j}}=0\ \forall j$ displays a superfluid to Mott-insulator transition for commensurate fillings $N/M$, and at a critical value of the ratio $U/t$ of interaction-to-tunnel energy. The phase boundaries of the transition are expected to be affected by the magnetic flux through an overall rescaling $t/U\to (t/U){\rm cos} (\Omega /M)$ [34]. The potential barrier considered here is localized on a single site j₀, i.e., ${{\Lambda }_{j}}=\Lambda {{\delta }_{j,{{j}_{0}}}}$ with ${{\delta }_{i,j}}$ being the Kronecker delta. As we will discuss in section 3.2, we find a superfluid-insulator transition even if the ring is interrupted by a weak link, although the phenomenon is rather a crossover, the ring being of finite size.

In this work, specific regimes of the system described by equation (1) will be captured analytically via the Tonks–Girardeau (TG) mapping (hard-core limit of infinite repulsions), and the mean-field Gross–Pitaevskii (GP) approximation (weak interactions and large fillings). To cover all the interaction regimes, numerical analysis will be also pursued, through truncated and exact diagonalization (ED) schemes and density-matrix renormalization-group (DMRG) methods. Details on these techniques are given in appendix C.

2.1. Identification of the qubit: effective two-level system

The Hamiltonian (1) is manifestly periodic in $\Omega$ with period $2\pi$. Therefore, we can restrict our study to the first rotational Brillouin zone, (actually to half of it, i.e., $\Omega \in [0,\pi ]$, due to the further symmetry $\Omega \leftrightarrow -\Omega$). In the absence of a barrier, the system is also rotationally invariant and therefore the particle–particle interaction energy does not depend on $\Omega$. The many-body ground-state energy, as a function of $\Omega$, is therefore given by a set of parabolas each corresponding to a well defined angular momentum state, shifted with respect to each other by a Galilean transformation and intersecting at the frustration points ${{\Omega }_{j}}=(2j+1)\pi$ [35, 36]. The presence of a finite barrier, $\Lambda \gt 0$, breaks the axial rotational symmetry and couples different angular momenta states, thus lifting the degeneracy at ${{\Omega }_{j}}$ by an amount $\Delta {{E}_{1}}$, see figure 1. The larger Λ, the larger is $\Delta {{E}_{1}}$, corresponding to the width of the gap separating the first two bands. Provided other excitations are energetically far enough from the two competing ground-states, this will identify the two-level system defining the desired qubit and its working point.

Zoom In Zoom Out Reset image size

Below, we discuss this issue with two different approaches: first, exploiting the mapping of the BH model to the quantum phase model, neglecting the fluctuations of the amplitude of the superfluid order parameter; this approach can capture, in particular, the regime of a large filling per lattice site [37, 38]. Then, via numerical calculation of the ground and first three excited energy levels of the BH model equation (1), we cover the case of lattice rings with a low filling.

Quantum phase model. In the regime of filling much larger than one, the number fluctuations on each site can be neglected and the behavior of the system is governed by the quantum phase model [37, 38] with Josephson couplings ${{J}_{j}}\sim \langle n\rangle |{{t}_{j}}|$, where $\langle n\rangle$ is the average number of bosons per well. The presence of a barrier constriction can be modeled by a weak link ${{J}_{{{j}_{0}}}}={{J}^{\prime }}\lt j$ ¹⁰ . The artificial flux $\Omega$ can be gauged away everywhere but at the j₀th site [39–41], thus giving rise to an energy term $-{{J}^{\prime }}{\rm cos} ({{\phi }_{{{j}_{0}}}}-{{\phi }_{{{j}_{0}}+1}}-\Omega )$. In this situation an effective action can be derived, which depends on a single phase difference $\theta ={{\phi }_{{{j}_{0}}}}-{{\phi }_{{{j}_{0}}+1}}$ across the weak link [25, 42]. The corresponding effective potential reads ${{V}_{{\rm eff}}}(\theta )=J{{\theta }^{2}}/M-J^{\prime} {\rm cos} (\theta -\Omega )$,which, for large $MJ^{\prime} /J$ and moderate M, defines a two-level system with degeneracy point at $\Omega =\pi$, as pictorially illustrated in figure 1. In the co-rotating frame, these two states correspond to the symmetric and antisymmetric combination of counter-circulating currents, where the degeneracy is split because of the inter-well tunneling.

BH model. Here we study the low-lying spectrum of the BH model (1) by a numerical analysis, performed in the dilute limit (low filling regime). This is complementary to the quantum phase model, in that we take into account the effect of the number fluctuations, and hence of the amplitude of the superfluid order parameter, on the lattice sites. In figure 2 we show the ED results for M = 16 and N = 4. The top-left panel shows how large interactions and moderate barrier strengths cooperate to define a doublet of energy levels at $\Omega =\pi$, well separated in energy with respect to the higher excited states; weaker interactions and larger barrier strengths, in contrast, do not allow for a clear definition of a two-level system (top-right panel). We observe that for increasing Λ, as expected, the gap increases and the bands become flatter, thus weakening the dependence of the energy on $\Omega$. The lower two panels display a complete analysis of the behavior of the spectral gap and its distance to the next excited level at $\Omega =\pi$ as a function of interactions and barrier strength, allowing us to identify the parameter regime for the existence of an effective two-level system. We notice in particular that weakly interacting gases cannot give rise to a sensible qubit within this approach, since one cannot isolate two levels out of the many-body spectrum with the sole tuning of the barrier strength, while this is possible for larger interaction strengths U.

Zoom In Zoom Out Reset image size

Using the above results, we conclude that the low-energy spectrum of the system (1) may define a qubit over a broad range of lattice filling values. It is vital for the manipulation of the qubit, though, to explore its quality. This implies in particular to study the dependence of $\Delta {{E}_{1}}$ on the system size and on the interaction strength, as will be considered in the next section 3. We will also analyze the nature of the qubit states; this will be the subject of section 4.

2.2. Density profiles

Before presenting our results concerning the quality of the qubit, we first focus on the density profiles of the gas close to the qubit working point, $\Omega =\pi$.

An evident effect of the barrier is a suppression of the particle density in its immediate proximity; depending on the ring size, the whole density profile along the ring may well be affected. The interplay between the interaction strength U and the barrier intensity $\Lambda$ implies different behaviors [43], as exemplified in figure 3 (panels (a)–(d)) for relatively small rings. The depth of the density depression increases monotonously with $\Lambda$ (inside each panel), while its width decreases with increasing U (see the different panels) since the density can be suppressed at the impurity site at the expense of multi occupancy of the other sites; the latter effect implies a non trivial dependence of the healing length on interaction strength. At strong repulsive interactions we also observe small Friedel-like oscillations of the density, which are a consequence of the peculiar strong correlations of 1D bosons that make their response to impurities similar to fermions.

Zoom In Zoom Out Reset image size

We note that, a sufficiently large barrier (at fixed U) makes the density profile vanish, thus effectively disconnecting the ring (panels (a)–(d) of figure 3). The barrier strength required to disconnect the ring depends on the interaction strength. Panel (e) of figure 3 shows the result of a thorough analysis of the transition line in the Λ-U plane: for a wide range of interaction strengths, the critical barrier height ${{\Lambda }_{{\rm c}}}$ displays a nearly perfect linear behavior with U. The prefactor turns out to be nearly proportional to the filling.

In this section we study in detail the spectroscopy of the qubit. We will analyze how the energy gaps $\Delta {{E}_{1}},\Delta {{E}_{2}}$ between the ground and, respectively, the first-excited / second-excited energy levels of the many-body Hamiltonian (1) depend on the system size and on the filling, for different $\Lambda$ and U. We find that the qubit is well resolved in the mesoscopic regime of intermediate ring sizes, and that it is at best separated from the higher energy levels of the many-body spectrum in the regime of strong interactions and weak barrier.

3.1. Scaling with the system size

In figure 4 we show both the qubit gap $\Delta {{E}_{1}}$ and the separation of the two levels from the rest of the spectrum in terms of $\Delta {{E}_{1}}/\Delta {{E}_{2}}$, as obtained by DMRG simulations at constant filling $N/L=1/4$ (see appendix C.4). The three panels correspond to different barrier intensities, from very weak to very high; each panel containing the three curves at varying interactions from moderate to hard-core. A clear power-law decay of $\Delta {{E}_{1}}$ results in all the regimes; the exponents depend on the interplay between the barrier and interaction strengths.

Zoom In Zoom Out Reset image size

In the small-barrier limit, we can work out the observed scaling law of the gap analytically resorting to the Luttinger-liquid effective field theory (see, e.g., [43]). Indeed we obtain that the quantum fluctuations of the density renormalize the barrier strength according to ${{\Lambda }_{{\rm eff}}}\sim \Lambda {{(d/L)}^{K}}$, where d is a short distance cut-off of the low-energy theory, $L=aM$ is the system size, a being the lattice spacing, and K is the Luttinger parameter [43]. This yields the scaling of the gap with M as

$$\begin{eqnarray}&&\Delta {{E}_{1}}\sim \nu {{\Lambda }_{{\rm eff}}}\sim {{M}^{-K}}\ ,\end{eqnarray} \tag{ 2 }$$

in agreement with the result found in [44] for a single impurity potential. As illustrated in panel (a) of figure 4, we find a very good agreement between the numerical data and the power-law behavior dictated by the Luttinger parameter obtained via the Bethe-Ansatz solution of the continuous model (a Lieb-Liniger gas [45]), suitable in the dilute limit of the BH model ¹¹ . For stronger barriers, interestingly, we observe in figures 4(b)–(c) that the gap scales again as a power-law, beyond the regime of validity of the analytical predictions. We also notice that the scaling of the gap is closely related to the scaling of the persistent currents flowing along the ring [47], which is determined by the shape of the ground state energy band.

By looking at the separation of the effective two-level system from the rest of the spectrum (dashed lines in figure 4), we can then start to identify an ideal regime of size, interaction and barrier for a realistic operational realization of the qubit. At low barrier intensity $\Lambda /t=0.1$ (panel (a)), indeed, a mesoscopic lattice of few tens of sites filled with mildly interacting bosons appears to be the best choice, since it would allow for a qubit gap of some ${{10}^{-3}}t$, while this being only a $\simeq {{10}^{-2}}$ fraction of the second excitation energy. Rings that are too large in size would improve the definition of the two-level system, yet at the price of too small a resolution of the qubit levels for practical addressing. When the barrier becomes stronger, the size dependence of $\Delta {{E}_{1}}/\Delta {{E}_{2}}$ becomes less and less important, with its absolute value increasing more and more (i.e., the qubit gets less and less isolated). Still, at intermediate barrier strengths $\Lambda /t=1$ (panel (b)), a nicely addressable pair of levels with splitting of some ${{10}^{-2}}t$, and a relative separation from the spectrum of order ${{10}^{-1}}t$, can be obtained in a mesoscopic lattice of $M\simeq 16$ sites with relatively weak interactions U = t. Conversely, if the barrier is strong enough to effectively cut the ring, the low lying levels of the many-body spectrum get almost equally spaced and therefore the qubit definition results to be poor.

3.2. Dependence of the qubit energy spectrum on the filling factor in mesoscopic rings

We concentrate next on the mesoscopic regime of few lattice sites, where, according to our scaling analysis at fixed small filling, the qubit enjoys simultaneously a clear definition with respect to the other excited states and a good energy resolution.

In figures 5 and 6 we present our results for the gap $\Delta {{E}_{1}}$ and $\Delta {{E}_{2}}$ as a function of the filling at fixed system size, studying its dependence on the barrier and on the interaction strength. The top panels of figure 5 present the data for fixed interaction strength ($U/t=1$ and $U/t=10$, respectively) with the curves representing barrier strengths from weak to strong. At small U (top-left), we observe a smooth dependence of $\Delta {{E}_{1}}$ on the boson filling, as expected in the superfluid regime of the Hamiltonian (1), of which the small ring is reminiscent. The increase of the barrier strength has two effects: first, at fixed filling, it increases the gap since it enhances the effect of the breaking of the rotational invariance and therefore lifts the degeneracy at half-flux. In addition, it changes the dependence of the gap on the filling from being monotonically decreasing to monotonically increasing, passing through a crossover situation. Since the healing length scales as $\xi \propto 1/\sqrt{\nu U}$, at small barrier strengths, a weakly interacting Bose gas screens the barrier, effectively reducing its height as the density is increased. On the other hand, for a large barrier, the system is effectively in the tunnel limit, and the situation is reversed. The barrier strength is effectively enhanced, since the tunnel energy required to move one particle from one side of the barrier to the other increases if the number of particles or the interaction strength are increased (in order to accommodate the tunnelling particle, the other particles have to readjust their configuration).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At large U (top-right) in figure 5, $\Delta {{E}_{1}}$ displays a more complex dependence on the filling, with pronounced peaks at particle numbers commensurate (or quasi) with the size, related to the presence of Mott lobes in the phase diagram of Hamiltonian (1) [48]. For weak barrier, indeed, the peaks appear at integer values of $N/M$, while for very strong potential constrictions the density is suppressed on one site: the system is close to a lattice with $M-1$ sites and peaks are consequently shifted. At intermediate barrier strengths we can observe a transient between the two regimes and broader peaks appear. Considering the very small system size, this effect arises because the presence of the healing length affects the whole bosonic density profile of the ring.

The top panels of figure 6 present data for fixed barrier strength ($\Lambda =0.1$ and $\Lambda =1$, respectively) with the curves representing interaction strengths from extremely weak to infinite values. First, we can clearly see the non-monotonous dependence of the gap on U, which was illustrated in figure 2, to hold at all fillings in both panels. Secondly, we notice that the dependence of $\Delta {{E}_{1}}$ on N drastically changes increasing the interaction strength, displaying different regimes: quickly decreasing, non-monotonous and almost constant. The rapid decrease of the energy gap at weak interactions can be understood (through a perturbative argument) in terms of level mixing of single-particle energies, which increases with the number of bosons involved [49]. In the opposite regime of hard core bosons, the energy gap is of the same order as the one of the non-interacting Fermi gas. This can be readily understood in terms of the TG Bose–Fermi mapping: indeed, in a non-interacting Fermi gas the energy gap is given by $\Delta {{E}_{1}}=(\sum _{j=1}^{N-1}{{\epsilon }_{j}}+{{\epsilon }_{N+1}})-\sum _{j=1}^{N}{{\epsilon }_{j}}$, where ${{\epsilon }_{j}}$ are the single-particle energies. In particular, for a small barrier, using perturbation theory, one obtains that the single-particle energy gaps ${{\epsilon }_{j+1}}-{{\epsilon }_{j}}$ are identical for all the avoided levels crossings, hence the gap $\Delta {{E}_{1}}$ is independent of the filling.

The lower panels of figures 5 and 6 display the ratio $\Delta {{E}_{1}}/\Delta {{E}_{2}}$. This allows us to identify the low-barrier, intermediate-to-large interaction regime at arbitrary filling as the most favorable for the qubit. Indeed, depending on the interaction strength, a too large barrier yields an unfavorable situation similar to the one depicted in the top-right panel of figure 2, where $\Delta {{E}_{2}}\sim \Delta {{E}_{1}}$. It is interesting to notice that these unfavorable cases correspond to values of barrier and interaction strength in the right panel of figure 3 where the ring is effectively disconnected. This allows us to identify the ratio $\Lambda /U$ as a useful parameter to define the quality of the qubit in terms of its energy resolution: the most advantageous parameter regime for the qubit corresponds to the lower half-plane in figure 3(e), below the critical line.

In summary, this analysis shows that a particularly favorable regime for the energy resolution of the qubit is the TG and small-barrier limit, where the system has a well defined gap, independent of the particle number and well separated from the remaining part of the many-body spectrum. However, for the realization of a tunable-gap qubit, the limits of weak interaction with low filling and intermediate interaction with high filling can be useful.

We close this section providing the order of magnitude for the gaps discussed above. For a ⁸⁷Rb gasin a mesoscopic ring shaped deep optical lattice of $\sim 50$ μm circumference and ten lattice wells, the hopping energy is of the order of $t\sim 0.5$ kHz. This yields a typical energy scale for the gap of tens to few 100 of Hz, depending on the choice of barrier strength, well within the range of experimental accessibility.

So far we focused our analysis on the behavior of the energy spectrum of the qubit as a function of the system parameters. We now investigate the ground state of the system in more detail. Special care is devoted to the regimes corresponding to a macroscopic superposition of circulation states. We assess the detectability of the latter through the study of the momentum distribution.

The momentum distribution is experimentally accessible in cold atoms experiments via TOF expansion measurements, by averaging over many repeated TOF realizations [50, 51], and is employed to get information about the current circulation along the ring [52–54]. It is defined as the Fourier transform with respect to the relative coordinate of the one-body density matrix ${{\rho }_{(1)}}({\bf x},{{{\bf x}}^{\prime }})=\langle {{\hat{\psi }}^{\dagger }}({\bf x})\hat{\psi }({{{\bf x}}^{\prime }})\rangle$:

$$\begin{eqnarray}&&n({\bf k})=\int {\rm d}{\bf x}\int {\rm d}{{{\bf x}}^{\prime }}\langle {{\hat{\psi }}^{\dagger }}({\bf x})\hat{\psi }({{{\bf x}}^{\prime }})\rangle {{{\rm e}}^{{\rm i}{\bf k}\cdot ({\bf x}-{{{\bf x}}^{\prime }})}},\end{eqnarray} \tag{ 3 }$$

where ${\bf x}$ and ${\bf x}^{\prime}$ denote the position of two points along the ring's circumference. Although, in general, ${\bf k}$ is a three-dimensional wave vector, here we restrict to consider a TOF picture along the symmetry axis of the ring, and therefore two-dimensional ${\bf k}$'s. To adapt equation (3) to our lattice system, we use $\hat{\psi }({\bf x})=\sum _{j=1}^{M}{{w}_{j}}({\bf x}){{\hat{b}}_{j}},$ where ${{w}_{j}}({\bf x})=w({\bf x}-{{{\bf x}}_{j}})$ is the Wannier function localized on the jth lattice site, and ${{{\bf x}}_{j}}$ denotes the position of the jth lattice site. Thereby, equation (3) can be recast into

$$\begin{eqnarray}&&n({\bf k})=|\tilde{w}({\bf k}){{|}^{2}}\mathop{\sum }\limits_{l,j=1}^{M}{{{\rm e}}^{{\rm i}{\bf k}\cdot ({{{\bf x}}_{l}}-{{{\bf x}}_{j}})}}\langle \hat{b}_{l}^{\dagger }{{\hat{b}}_{j}}\rangle ,\end{eqnarray} \tag{ 4 }$$

where $\tilde{w}({\bf k})$ is the Fourier transform of the Wannier function.

To avoid effects of the proximity of the superfluid-insulator transition, in the following analysis, we focus on incommensurate fillings (see section 3.2 for a more detailed discussion).

In absence of barrier $\Lambda =0$, the system has no circulation for $\Omega \lt \pi$ and one quantum of circulation for $\Omega \gt \pi$, while at the frustration point $\Omega =\pi$, it is a perfectly balanced superposition of the two states. As a consequence, the momentum distribution is peaked at ${\bf k}=0$ for $\Omega \lt \pi$ and is ring-shaped for $\Omega \gt \pi$, as discussed in appendix A. At $\Omega =\pi$, instead, it displays an interference of the two situations, reflecting the coherent superposition of the two states (see appendix B). When $\Lambda \ne 0$, the superposition state occurs for a wide range of $\Omega$, thereby displaying interference effects as shown in figure 7. The relative weight of the two-quanta-of-circulation components in the superposition strongly depends on $\Omega$, $\Lambda$, and U. In particular, at the frustration point, the superposition is perfectly balanced, independently of $\Lambda$, and U. Away from the frustration points, the relative weights tend to the unperturbed ones carrying zero or one quantum of angular momentum. This phenomenon occurs over a distance in $\Omega -\pi$ that depends on $\Lambda$: the smaller is $\Lambda$, the faster the unperturbed weights are recovered. For this reason, in figures 7 and 9, we slightly off-set $\Omega$ from the frustration point (the weights of the circulating states are not equal, yet close enough to ensure that both angular momentum states contribute significantly to the superposition). For $\Omega \gt \pi$, the component carrying one quantum of angular momentum has a larger weight in the superposition, making the effect of the barrier and its screening easily detectable in the TOF image; the opposite situation occurs for $\Omega \lt \pi$. The TOF results shown in figures 7 and 9 quantitatively depend on the choice of $\Omega$, but the screening effect of the barrier and the detectability of the superposition result weakly affected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To understand the TOF results of figure 7, it is instructive to consider first the case without interactions, U = 0, that is analytically accessible. The corresponding momentum distribution close the frustration point and for a weak barrier reads (see equation (B.4) in appendix B for the derivation)

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} n({\bf k}) & = & {{{\rm sin} }^{2}}(\varphi /2)J_{n}^{2}(|{\bf k}|R)+{{{\rm cos} }^{2}}(\varphi /2)J_{n+1}^{2}(|{\bf k}|R) \\ {} & {} & +{\rm sin} (\varphi ){\rm cos} ({{\gamma }_{{\bf k}}}){{J}_{n}}(|{\bf k}|R){{J}_{n+1}}(|{\bf k}|R), \\ \end{array}\end{eqnarray} \tag{ 5 }$$

where the Bessel functions J_n correspond to states with angular momentum n, and ${{\gamma }_{{\bf k}}}$ is the angle along the ring; the parameter φ is a function of the flux and the barrier strength (see equation (B.3) in appendix B). Equation (5) shows that the TOF images allow to visualize the superposition between of states with different angular momenta: the functions J₀ and J₁ interfere, giving rise to a peak at zero ${\bf k}$ and a fringe with ring-shaped symmetry. The detectability of this feature increases with the barrier strength $\Lambda$. Note that the angular position of the peak in momentum space depends on the position of the barrier in real space along the ring; it would be affected by a phase shift between the two states of well-defined angular momentum.

The superposition state for small U can be analyzed in a similar way. We note in figure 7 that, for sufficiently weak interactions, an angular modulation of the ring-shaped momentum distribution arises. A stronger barrier makes the angular asymmetry increasing, while the interaction strength, by screening the barrier, leads to the opposite phenomenon.

Upon increasing the interaction strength from intermediate to very large, we observe a smearing of the modulated ring shape TOF images. This is an effect of increased quantum fluctuations, which leads, for strong barrier strengths, to a single maximum centred at non-zero ${\bf k}$ values.

The very different TOF images between the regimes of weak and strong interactions can be understood by recalling the different nature of the superposition state in the various interaction regimes [27, 29]. For instance, at zero or very weak interactions, within the GP regime, the many-body state is a coherent state of single particle superpositions. Increasing the interaction strength to the intermediate regime the superposition is described by the so-called NOON state $|N,0\rangle +|0,N\rangle$, i.e., a macroscopic superposition of states where all bosons occupy either the state with zero angular momentum or the one carrying one quantum of angular momentum. For increasing interactions this many-body entangled state matches the known macroscopic superposition of Fermi spheres at very large interactions [27, 55].

For all regimes of interactions, we notice that the TOF images become independent of the barrier above a critical value of the barrier strength, which well agrees with the critical value ${{\Lambda }_{{\rm c}}}$ for disconnecting the ring, as identified in figure 3. Globally, we observe that good-quality TOF images allowing to easily identify the superposition of current states as a modulated ring structure are found for a ratio $\Lambda /U$ in the vicinity or above the critical line of figure 3(e).

To quantify the detectability of the superposition state in the momentum distribution for different barrier and interaction strengths, we define the averaged contrast between the momentum distribution with and without barrier

$$\begin{eqnarray}&&\eta =\frac{\int {\rm d}{\bf k}|{{n}_{\Lambda \ne 0}}({\bf k})-{{n}_{\Lambda =0}}({\bf k})|}{\int {\rm d}{\bf k}\ {{n}_{\Lambda \ne 0}}({\bf k})+{{n}_{\Lambda =0}}({\bf k})},\end{eqnarray} \tag{ 6 }$$

reflecting the modification in the integrated momentum distribution due to superposition of states induced by the barrier. We find that $\eta$ is non-monotonic upon increasing the interactions between the particles, while keeping fixed the barrier strength—figure 8. This is an effect of the non-monotonic screening of the barrier as a function of interaction strength, first predicted in [43] through the study of the persistent-current amplitude.

Finally, we comment on the expected behavior for a system with filling larger than one. In figure 9 we show the TOF images for larger fillings, ranging from values of $N/M$ close to one, obtained with truncated ED, to fillings much larger than one, obtained solving the GP equation (C.3). In both cases we note that the TOF images are qualitatively the same as the ones shown in figure 7, and therefore our analysis is relevant also for systems with larger number of particles, like the ones employed in the experiments so far. We notice that, at higher filling, a larger barrier strength is needed, with respect to the lower filling case, to produce the same superposition and to observe the same TOF.

We considered a system of bosonic atoms loaded in a ring-shaped 1D optical lattice potential, hosting a localized barrier on a given site of the lattice. A ring lattice could be produced, for example, through protocols based on interference patterns of Laguerre–Gauss or Bessel laser beams [7, 14] or through spatial light modulators [25]. Besides its possible exploitation in quantum technology, the system provides a paradigmatic arena to study the interplay between quantum fluctuations, interactions and the role of the barrier potential.

In order to address the qubit effective quantum dynamics encoded into the system, its coherent superposition of current states and the scaling of its properties with system size, we considered ring sizes ranging from few lattice sites ($10-20$) to larger structures (100). We provide a direct evidence that the qubit dynamics can be achieved beyond the pure superfluid phase dynamics conditions (described by the quantum phase model) exploited to derive the effective double-well potential [25], see figures 1 and 2. Incidentally, we note that one- and two-qubit gates can be realized with our quantum device by tuning the the barrier height and interaction suitably (see [25] for details).

We studied the resolution and detectability of the qubit, by following two routes: the analysis of the scaling properties (both with the number of particles and system size) of the energy gap between the two energy states, and the analysis of the momentum distribution.

The energy gap of the two-level system. We quantified the scaling of the gaps $\Delta {{E}_{1}}$ and $\Delta {{E}_{2}}$ with system size. Our results indicates that $\Delta {{E}_{1}}$ is appreciable for small and mesoscopic systems while it is suppressed in thermodynamic limit (figure 4), decaying as a power law with system size. This follows from the localization of the barrier to a point like on the scale of the lattice spacing, and shows that the lattice potential along the ring bears several added values with respect to the uniform ring case or the case of a broad barrier [27], ultimately facilitating the exploitation of the device in future atomtronic integrated circuits. Our scaling analysis allows us to identify the mesoscopic regime as the most suited for the realization of the AQUID—see figure 4.

Momentum distribution. For the mesoscopic structures, we demonstrated that the coherent superposition of forward and backward scattering of the particles through the barrier site, is indeed detectable through TOF expansion. This is a strong a posteriori evidence of the two-level-system effective physics that is encoded into the system. For fixed values of the filling parameter, the detectability of the superposition depends on the relative size between barrier and interaction strengths: the barrier makes the detectability increasing, while the interaction strength, screening the barrier, leads to the opposite phenomenon, yielding a non-monotonous behavior of the detectability and averaged contrast. By increasing the filling parameter for fixed U and $\Lambda$, the screening of the barrier is enhanced, and therefore the barrier is less effective in creating the coherent superposition of flows. A separate discussion is needed for the regime of large interactions. This regime is characterized by fermionic effects due to strong correlations. In particular, we note a good detectability of the superposition state with a simultaneous presence of density modulations along the ring.

In summary, our work indicates that a mesoscopic ring lattice with localized barrier provides a candidate for the AQUID. The qubit dynamics is detectable in a wide range of system parameters. We have identified the ratio $U/\Lambda$ as an important parameter to discuss the behavior of the qubit both in terms of its definition with respect to the rest of the many-body spectrum (the qubit turns out to be best defined below the critical line in figure 3(e)) and its detectability in TOF images (the contrast is found to be best defined around or above the same critical line). This allows us to conclude that the ratio $U/\Lambda$ on the critical line yields an optimal parameter choice, corresponding to the best trade-off between the momentum distribution detectability and the gap resolution. Given the flexibility achieved in the actual experiments, we believe that such an optimum parameter regime is well within the current experimental know-how of the field.

DA is indebted to S Vinjanampathy and PN Ma for help implementing the exact diagonalization algorithm. We thank F Auksztol, H Crepaz, and R Dumke for discussions. We acknowledge support from the Merlion project 'LUMATOM', the Institut Universitaire de France, the ERC Handy-Q grant N.258608, the ANR project Mathostaq ANR-13-JS01-0005-01, and the Italian FIRB project RBFR12NLNA. DMRG numerical simulations have been performed on the MOGON cluster of JGU-Mainz.

The signature of a non-vanishing current flow along the ring lattice is a ring-shaped configuration of the momentum distribution (see appendix B for a derivation in the non-interacting limit). Figure A1 shows the predicted TOF images in the absence of the barrier for various interaction strengths. The perfect ring shape reflects angular momentum conservation at all interaction strengths, consistent with Leggett's theorem establishing that the persistent currents through a rotationally invariant system are not affected by the interactions. We note, however, that the detectability in the TOF images is reduced at large interactions, due to the enhanced role of phase fluctuations.

Zoom In Zoom Out Reset image size

In the non-interacting regime the many-body problem reduces to a single-particle one. In the absence of the barrier the Schrödinger equation, in polar coordinates and scaling the energies in units of ${{E}_{0}}={{\hbar }^{2}}/2\;m{{L}^{2}}$, with m being the particle mass, and L the system size, reads

$$\begin{eqnarray*}&&{{\left( -{\rm i}\frac{\partial }{\partial \theta }-\frac{\Omega }{2\pi } \right)}^{2}}\psi (\theta )=E\psi (\theta ),\end{eqnarray*}$$

where $\theta \in [0,2\pi ]$. The wavefunction for a state with defined angular momentum is a plane wave $\psi (\theta )=(1/\sqrt{2\pi }){{{\rm e}}^{{\rm i}n\theta }}$, where $n\in \mathbb{Z}$ to satisfy periodic boundary conditions, and the corresponding spectrum is ${{E}_{n}}={{(n-\Omega /2\pi )}^{2}}$. The momentum distribution then reads

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} n({\bf k}) & = & \int {\rm d}{\bf x}\int {\rm d}{\bf x}^{\prime} {{{\rm e}}^{{\rm i}{\bf k}\cdot ({\bf x}-{{{\bf x}}^{\prime }})}}{{\psi }^{*}}({\bf x})\psi ({{{\bf x}}^{\prime }}) \\ {} & \sim & {{\left| \int _{0}^{2\pi }{\rm d}\theta \ {{{\rm e}}^{{\rm i}({{k}_{x}}R{\rm cos} \theta +{{k}_{y}}R{\rm sin} \theta )}}{{\psi }^{*}}(\theta ) \right|}^{2}} \\ {} & = & {{\left| {{{\rm e}}^{{\rm i}m\gamma }}{{J}_{m}}(|{\bf k}|R) \right|}^{2}}={{\left| {{J}_{n}}(|{\bf k}|R) \right|}^{2}}, \\ \end{array}\end{eqnarray} \tag{ B.1 }$$

where $R=L/2\pi$ is the ring radius, we have defined γ as ${{k}_{x}}=|{\bf k}|{\rm sin} \gamma$, ${{k}_{y}}=|{\bf k}|{\rm cos} \gamma$, and J_n is the nth order Bessel function of the first kind. For n = 0 the momentum distribution is peaked at ${\bf k}=0$, while for $n\gt 0$ it is ring shaped, with a radius that grows with n.

In the presence of a localized barrier of strength λ the Schrödinger equation becomes

$$\begin{eqnarray*}&&{{\left( -{\rm i}\frac{\partial }{\partial \theta }-\frac{\Omega }{2\pi } \right)}^{2}}\psi (\theta )+\lambda \delta (\theta )\psi (\theta )=E\psi (\theta ).\end{eqnarray*}$$

The effect of the δ-barrier is to mix states with different angular momentum. For a small barrier we can reduce to the simplest case of mixing of states that differ by just one quantum of angular momentum, and apply degenerate perturbation theory. We write the Hamiltonian in the following form

$$\begin{eqnarray}H=\left( \begin{array}{ccccccccccccccc} {{E}_{n}} & \lambda /2\pi \\ \lambda /2\pi & {{E}_{n+1}} \\ \end{array} \right);\end{eqnarray} \tag{ B.2 }$$

the corresponding eigenvalues and eigenvectors reads:

$$\begin{eqnarray*}&&{{\epsilon }_{1,2}}=\frac{{{E}_{n+1}}+{{E}_{n}}}{2}\pm \frac{\sqrt{\delta {{E}^{2}}+{{\lambda }^{2}}/{{\pi }^{2}}}}{2}\ ,\end{eqnarray*}$$

where $\delta E={{E}_{n+1}}-{{E}_{n}}$, and

$$\begin{eqnarray*}&&{{w}_{1}}=\left( \begin{array}{ccccccccccccccc} {\rm sin} (\varphi /2) \\ {\rm cos} (\varphi /2) \\ \end{array} \right)\;,\ \ \ \ \ \ \ \ {{w}_{2}}=\left( \begin{array}{ccccccccccccccc} {\rm cos} (\varphi /2) \\ -{\rm sin} (\varphi /2) \\ \end{array} \right)\ ,\end{eqnarray*}$$

where

$$\begin{eqnarray}&&{{{\rm cos} }^{2}}(\varphi /2)=\frac{\sqrt{\delta {{E}^{2}}+{{\lambda }^{2}}/{{\pi }^{2}}}-\delta E}{2\sqrt{\delta {{E}^{2}}+{{\lambda }^{2}}/{{\pi }^{2}}}}\ .\end{eqnarray} \tag{ B.3 }$$

We then write the wavefunction as

$$\begin{eqnarray*}&&\psi (\theta )=\frac{1}{\sqrt{2\pi }}{\rm sin} (\varphi /2){{{\rm e}}^{{\rm i}n\theta }}+\frac{1}{\sqrt{2\pi }}{\rm cos} (\varphi /2){{{\rm e}}^{{\rm i}(n+1)\theta }},\end{eqnarray*}$$

where $\varphi$ depends on $\lambda$ and $\Omega$.

The momentum distribution in this case becomes

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} n({\bf k}) & \sim & {{\left| \int _{0}^{2\pi }{\rm d}\theta \ {{{\rm e}}^{{\rm i}({{k}_{x}}R{\rm cos} \theta +{{k}_{y}}R{\rm sin} \theta )}}{{\psi }^{*}}(\theta ) \right|}^{2}} \\ {} & = & {{\left| {\rm sin} (\varphi /2){{{\rm e}}^{{\rm i}n\gamma }}{{J}_{n}}(|{\bf k}|R)+{\rm cos} (\varphi /2){{{\rm e}}^{{\rm i}(n+1)\gamma }}{{J}_{n+1}}(|{\bf k}|R) \right|}^{2}} \\ {} & = & {{{\rm sin} }^{2}}(\varphi /2)J_{n}^{2}(|{\bf k}|R)+{{{\rm cos} }^{2}}(\varphi /2)J_{n+1}^{2}(|{\bf k}|R) \\ {} & {} & +2{\rm sin} (\varphi /2){\rm cos} (\varphi /2){\rm cos} (\gamma ){{J}_{n}}(|{\bf k}|R){{J}_{n+1}}(|{\bf k}|R)\ , \\ \end{array}\end{eqnarray} \tag{ B.4 }$$

where an interference term, proportional to ${\rm cos} \gamma$, appears between the two states with defined angular momentum, giving rise to a $2\pi$-periodic angular modulation of the ring shape found previously. This behavior is the same found in figure 7, where we observe an analogous modulation in the weak barrier and weak interaction case, that we can interpret than as direct consequence of the superposition of two stated that differ by one quantum of angular momentum.

C.1. Infinite interaction limit: TG gas

In the limiting case of infinite repulsive contact interaction between the particles ($U\to \infty$), the so-called hard-core bosons or TG gas, an exact approach can be pursed to diagonalize Hamiltonian (1). Since multi-occupancy of one site is forbidden by the infinite interaction energy, it can be simplified into

$$\begin{eqnarray}&&H=-t\mathop{\sum }\limits_{i=1}^{M}({{{\rm e}}^{-{\rm i}\Omega /M}}b_{i}^{\dagger }{{b}_{i+1}}+{\rm h}.{\rm c}.)+t\mathop{\sum }\limits_{i=1}^{M}{{\Lambda }_{i}}n_{i}^{b},\end{eqnarray} \tag{ C.1 }$$

where the bosonic annihilation and creation operators have the additional on-site constraints $b_{i}^{2}=b_{i}^{\dagger 2}=0$ and ${{b}_{i}}b_{i}^{\dagger }+b_{i}^{\dagger }{{b}_{i}}=1$. By applying the Jordan–Wigner transformation

$$\begin{eqnarray*}&&{{b}_{j}}=\Pi _{l=1}^{j-1}{{{\rm e}}^{{\rm i}\pi f_{l}^{\dagger }{{f}_{l}}}}{{f}_{j}}\ ,\end{eqnarray*}$$

where f_i ($f_{i}^{\dagger }$) are fermionic annihilation (creation) operators, the Hamiltonian (C.1) can be mapped into the one for spineless fermions:

$$\begin{eqnarray}&&H=-t\mathop{\sum }\limits_{i=1}^{M}({{{\rm e}}^{-{\rm i}\Omega /M}}f_{i}^{\dagger }{{f}_{i+1}}+{\rm h}.{\rm c}.)+t\mathop{\sum }\limits_{i=1}^{M}{{\Lambda }_{i}}n_{i}^{f}\end{eqnarray} \tag{ C.2 }$$

This Bose–Fermi mapping is the analogous, for a discrete system, of the one introduced by Girardeau for a continuous system [56]. Hamiltonians (C.1) and (C.2) have the same spectrum, but non-trivial differences appear in the off-diagonal correlation functions: $\langle f_{i}^{\dagger }{{f}_{j}}\rangle$ versus $\langle b_{i}^{\dagger }{{b}_{j}}\rangle$, which we have calculated following the same procedure described in [57]. Such difference affects, for example, the momentum distribution, which is much narrower for hard-core bosonic systems than for non-interacting fermions. The density, and all the quantities related to it, are instead identical, see for example figure 3. This 1D peculiar strongly correlated TG phase has been demonstrated in several experiments on bosonic wires [58, 59].

C.2. GP equation

In the limit of weak interactions, we adopt a mean-field approximation to simplify the many-body Schrödinger equation. This is the GP equation for the bosons subjected to a lattice potential, in the presence of a gauge field:

$$\begin{eqnarray}&&\frac{{{\hbar }^{2}}}{2m}{{\left( -{\rm i}{{\partial }_{x}}-\frac{\Omega }{L} \right)}^{2}}\Psi +{{U}_{b}}\delta (x)\Psi +{{V}_{0}}{{{\rm sin} }^{2}}\left( \frac{M}{L}\pi x \right)\Psi +{{g}_{1{\rm D}}}|\Psi {{|}^{2}}\Psi =\mu \Psi ,\end{eqnarray} \tag{ C.3 }$$

where $\Psi$ is the condensate wavefunction, $\mu$ is the chemical potential, V₀ is the optical-lattice depth, U_b is the strength of the localized barrier, modeled as $\delta (x)$ in this continuous model, m is the particle mass, and g_1D is the effective interaction coupling strength in one-dimension, related to the three-dimensional scattering length a as ${{g}_{1{\rm D}}}=2{{\hbar }^{2}}a/ma_{\bot }^{2}$ [60].

The continuous-model barrier strength U_b is connected to the discrete-model one $\Lambda$ by ${{U}_{{\rm b}}}=\Lambda L/M$. In absence of the lattice potential, an analytical soliton solution for equation (C.3) has been found in [43]. In the further limiting case of vanishing interaction and small barrier strength, the expression for the wavefunction can be obtained perturbatively with respect to the barrier strength. This approach helps the understanding of the corresponding momentum distribution and TOF images (see appendix B). In the presence of the lattice potential, we solve equation (C.3) numerically by integrating it in imaginary times. We pursue this approach as a benchmark case for the BH model at weak interaction. Moreover the GP equation is a particularly suitable tool for the large-N regime, which is routinely realized in experiments.

C.3. ED schemes

C.3.1. Working in the full Hilbert space

The ED is a computational method in many-body physics [61, 62] which gives exact eigenstates and eigenvalues of the Hamiltonian without making any simplifying assumptions about the physical system. However the method is applicable to small systems and small fillings $N/M$. The reason for that is provided by the fact that the Hilbert space spanned by the many-particle Fock states cannot be too large. Specifically, to implement the ED, one has to consider all the possible combinations of N particles over M sites, modulo the permutations of identical particles. The dimension of the Hilbert space is given by [63]:

$$\begin{eqnarray}&&D=\frac{(M+N-1)!}{(M-1)!N!}\end{eqnarray} \tag{ C.4 }$$

In section 4, we considered values of the filling: $5/11,15/11,24/9$. Correspondingly, the Hilbert space dimensions for that fillings are $3003;3268760;10\;518\;300$.

The non-diagonal part of the BH model can be written efficiently with the help of sparse matrices routine. The ground and the first excited state state eigenvalues and eigenvectors of the system can be found explicitly with help of Lanczos algorithm [61, 62].

C.3.2. Working in the truncated Hilbert space

To study the system for larger sizes and larger fillings, the ED scheme works upon reducing the dimension of the Hilbert space in a controlled way. This is achieved by restricting (truncating) the lattice site occupation number up to some given integer number K. The main difficulty of the truncated ED is the generation of the truncated Hilbert space in a numerically efficient way ¹² . Here we are using the following algorithm to achieve the goal. At first we write the function $f(M,n,K)$ which splits a positive integer M into a sum of n positive integers (where each of the integers is smaller than K) up to commutativity (so this function is returning matrix). Then we define the number s, where $s=[M/K]+1$, if $M/K-[M/K]\gt 0$ and else $s=[M/K]$. Then the truncated Hilbert space can be generated in the following three steps: (1) to apply function $f(M,n,K)$ by changing n from s to N with a step 1; (2) to concatenate each line in the matrix which return $f(M,n,K)$ with required amount of zeros to make lines of matrix N-dimensional arrays; (3) to generate all possible permutations for any line of the matrix. The dimension of the truncated Hilbert space is given by the following expression [64]:

$$\begin{eqnarray}&&{{D}_{K}}=\mathop{\sum }\limits_{i=1}^{\left[ \frac{N}{K+1} \right]}{{(-1)}^{j}}\left( \begin{array}{ccccccccccccccc} M+N-1-j(K+1) \\ M-1 \\ \end{array} \right)\left( \begin{array}{ccccccccccccccc} M \\ j \\ \end{array} \right)\end{eqnarray} \tag{ C.5 }$$

where the brackets $[]$ stand for the floor function.

For example for the case of the filling $24/9$ and K = 6 which is considered in section 4, ${{D}_{6}}=2345\;553$ which is almost 4.5 times smaller the dimension of the full Hilbert space. Indeed, in this way one can reduce D for the several order of magnitudes, but that will introduce errors, especially at small U. Here we estimate the errors in the following way. We calculate the particle number fluctuations (variance) per lattice site

$$\begin{eqnarray}&&{{\sigma }_{i}}=\sqrt{\langle n_{i}^{2}\rangle -{{\langle {{n}_{i}}\rangle }^{2}}}\end{eqnarray} \tag{ C.6 }$$

We assume that: if $\langle {{n}_{i}}\rangle +5{{\sigma }_{i}}\lt K,\forall i$ then error (in calculating expectation values) is smaller than $0.0006\%$, if $\langle {{n}_{i}}\rangle +4{{\sigma }_{i}}\lt K$ then error is smaller than $0.006\%$, if $\langle {{n}_{i}}\rangle +3{{\sigma }_{i}}\lt K$ then error is smaller than 3% and if $\langle {{n}_{i}}\rangle +2{{\sigma }_{i}}\lt K$ for any i then error is smaller than 5%. For the filling $24/9$ an estimated error for U = 10 is less than $0.006\%$ and for U = 1 it becomes 5%.

The momentum distribution in figure 9 is calculated with the scheme detailed above.

C.4. DMRG method