Coherent superposition of current flows in an Atomtronic Quantum Interference Device

We consider a correlated Bose gas tightly confined into a ring shaped lattice, in the presence of an artificial gauge potential inducing a persistent current through it. A weak link painted on the ring acts as a source of coherent back-scattering for the propagating gas, interfering with the forward scattered current. This system defines an atomic counterpart of the rf-SQUID: the atomtronics quantum interference device (AQUID). The goal of the present study is to corroborate the emergence of an effective two-level system in such a setup and to assess its quality, in terms of its inner resolution and its separation from the rest of the many-body spectrum, across the different physical regimes. In order to achieve this aim, we examine the dependence of the qubit energy gap on the bosonic density, the interaction strength, and the barrier depth, and we show how the superposition between current states appears in the momentum distribution (time-of-flight) images. A mesoscopic ring lattice with intermediate-to-strong interactions and weak barrier depth is found to be a favorable candidate for setting up, manipulating and probing a qubit in the next generation of atomic experiments.


I. INTRODUCTION
The progress achieved in optical micro-fabrication has laid much the foundation of atomtronics: Bose-Einstein condensates manipulated in optical circuits of various spatial shapes, with lithographic precision [1][2][3][4][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][8][9][10][11][12][13][14]. 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 [15,16], 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 [17,18]. 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 a very focused blue-detuned laser creates [19][20][21][22] an effective potential constriction (giving rise to the aforementioned weak link).
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 [23,24]. 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). This issue was first considered for homogeneous rings, in the absence of a lattice: the qubit splitting, i.e., the spectral gap between the many-body first-excited and ground-state energies, has been shown to dramatically depend on the spatial form of the potential constriction. For a perfectly localized delta-barrier potential it was demonstrated that the gap scales down as a power of the ring length [25] at strong interactions. Conversely, a finite width of the barrier leads to an exponential suppression of the same gap with increasing bosonic density [26]. This is the situation for the above-mentioned experimental realizations [17][18][19][20][21][22].
We concentrate on the quality of the qubit in the cold-atom ring lattice: in particular, we characterize the energy structure around the original degeneracy point at half-flux quantum. By employing a combination of analytical and numerical techniques, we cover all the relevant physical regimes of interactions and fillings. We show that: i) the gap ∆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 ∆E 1 displaying a favorable dependence in a wide range of system parameters; iii) the superposition state is visible in the momentum distribution of the bosonic gas, 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 Sec. III the energy spectrum of the system and its scaling with system size, filling, and interaction is analyzed. In Sec. IV 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 Sec. V. Technicalities on the employed methods and further details are provided in the Appendices.

II. THE PHYSICAL SYSTEM
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 Ω, which can be experimentally induced for neutral atoms as a Coriolis flux by rotating the lattice at constant velocity [21,30], or as a synthetic gauge flux by imparting a geometric phase directly to the atoms via suitably designed laser fields [31][32][33].
In the tight-binding approximation, this system is described by the 1D Bose-Hubbard (BH) Hamiltonian where b j (b † j ) are bosonic annihilation (creation) operators on the jth site and n j = b † j b j is the corresponding number operator. Periodic boundaries are imposed by taking b M +1 ≡ b 1 . The parameter U takes into account the finite scattering length for the atomic two-body collisions on the same site; Λ 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 Ω is taken into account through the Peierls substitution: t → e −iΩ M (t is the hopping amplitude). In the thermodynamic limit, the BH model for Λ j = 0∀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 -π -π/2 0 π/2 π θ 0.8 Main panel: sketch of the qubit energy splitting, due to the barrier Λ, for the two lowest-lying energy states in the many-body spectrum of model (1). Black dashed lines denote the ground-state energy in absence of the barrier, as a function of the flux Ω. Switching on the barrier opens a gap at the frustration point Ω = π (continuous red lines). The three insets show the qualitative form of the effective potential at Ω = 0, π, 2π. Note the characteristic double-well shape forming at Ω = π. The qubit, or effective two-level system, corresponds to the two lowest energy levels of this potential. In this figure the energies are plotted in arbitrary units.
energy. The phase boundaries of the transition are expected to be affected by the magnetic flux through an overall rescaling t U → (t U ) cos(Ω M ) [34]. The potential barrier considered here is localized on a single site j 0 , i.e., Λ j = Λδ j,j0 with δ i,j being the Kronecker delta. As we will discuss in Sec. III B, 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 Eq. (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.
A. Identification of the qubit: effective two-level system The Hamiltonian (1) is manifestly periodic in Ω with period 2π. Therefore, we can restrict our study to the first rotational Brillouin zone, (actually to half of it, i.e., Ω ∈ [0, π], due to the further symmetry Ω ↔ −Ω). In the absence of a barrier, the system is also rotationally invariant and therefore the particle-particle interaction energy does not depend on Ω. The many-body ground-state energy, as a function of Ω, 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 inter-secting at the frustration points Ω j = (2j + 1)π [35,36]. The presence of a finite barrier, Λ > 0, breaks the axial rotational symmetry and couples different angular momenta states, thus lifting the degeneracy at Ω j by an amount ∆E 1 , see Fig. 1. The larger Λ, the larger is ∆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.
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 Eq. (1), we cover the case of lattice rings with a low filling.
a. 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 ∼ ⟨n⟩ t j , where ⟨n⟩ is the average number of bosons per well. The presence of a barrier constriction can be modeled by a weak link J j0 = J ′ < J [39] The artificial flux Ω can be gauged away everywhere but at the j 0 th site [40][41][42], thus giving rise to an energy term −J ′ cos(φ j0 − φ j0+1 − Ω). In this situation an effective action can be derived, which depends on a single phase difference θ = φ j0 − φ j0+1 across the weak link [24,43]. The corresponding effective potential reads V eff (θ) ≐ Jθ 2 M − J ′ cos(θ − Ω), which, for large M J ′ J and moderate M , defines a two-level system with degeneracy point at Ω = π, as pictorially illustrated in Fig. 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 interwell tunneling.
b. Bose-Hubbard 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 Fig. 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 Ω = π, 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 Ω. The lower two panels display a complete analysis of the behavior of the spectral gap and its distance to the next excited level at Ω = π 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 . 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 ∆E 1 on the system size and on the interaction strength, as will be considered in the next Sec. III. We will also analyze the nature of the qubit states; this will be the subject of Sec. IV.

B. 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, Ω = π.
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 Λ implies different behaviors [44], as exemplified in Fig. 3 (panels a-d) for relatively small rings. The depth of the density depression increases monotonously with Λ (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.  We note that, a sufficiently large barrier (at fixed U ) makes the density profile vanish, thus effectively disconnecting the ring (panels a-d of Fig. 3). The barrier strength required to disconnect the ring depends on the interaction strength. Panel e) of Fig. 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 Λ c displays a nearly perfect linear behavior with U . The prefactor turns out to be nearly proportional to the filling.

III. ENERGY GAP OF THE TWO-LEVEL CURRENT-FLOW SYSTEM
In this section we study in detail the spectroscopy of the qubit. We will analyze how the energy gaps ∆E 1 , ∆E 2 between the ground and, respectively, the first-excited / secondexcited energy levels of the many-body Hamiltonian (1) depend on the system size and on the filling, for different Λ and U . We shall see 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 manybody spectrum in the regime of strong interactions and weak barrier.
A. Scaling with the system size In Fig. 4 we show both the qubit gap ∆E 1 and the separation of the two levels from the rest of the spectrum in terms of ∆E 1 ∆E 2 , as obtained by DMRG simulations at constant filling N L = 1 4 (see App. 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 ∆E 1 results in all the regimes; the exponents depend on the interplay between the barrier and interaction strengths.
In the small-barrier limit, we can work out the observed scaling law of the gap analytically resorting to the Luttingerliquid effective field theory (see, e.g., Ref. [44]). Indeed we obtain that the quantum fluctuations of the density renormalize the barrier strength according to Λ eff ∼ Λ(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 [44]. This yields the scaling of the gap with M as, in agreement with the result found in Ref. [45] for a single impurity potential. As illustrated in panel a) of Fig. 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 [46]), suitable in the dilute limit of the BH model [47]. For stronger barriers, interestingly, we observe in Fig. 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 [48], 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 Fig. 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 Λ 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 ≃ 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 ∆E 1 ∆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 Λ 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 ≃ 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.
B. 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 Figs. 5 and 6 we present our results for the gap ∆E 1 and ∆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 Fig. 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 ∆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 stronger breaks the rotational invariance and therefore lifts the degeneracy at halfflux. In addition, it changes the dependence of the gap on the filling from being monotonically decreasing to monotonically increasing, passing through a crossover situation. At small barrier strengths, the gap decreases due to a barrier renormalization with filling. At large U (top-right) in Fig. 5, ∆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) [49]. 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 Fig. 6 present data for fixed barrier strength (Λ = 0.1 and Λ = 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 Fig. 2, to hold at all fillings in both panels. Secondly, we notice that the dependence of ∆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 [50]. In the opposite regime of hard core bosons, the energy gap is of the same order as the one of the non-interacting Bose 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 where j are the singleparticle energies. In particular, for a small barrier, using perturbation theory, one obtains that the single-particle energy gaps j+1 − j are identical for all the avoided levels crossings, hence the gap ∆E 1 is independent of the filling.
The lower panels of Figs. 5 and 6 display the ratio ∆E 1 ∆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 Fig. 2, where ∆E 2 ∼ ∆E 1 . It is interesting to notice that these unfavorable cases correspond to values of barrier and interaction strength in the right panel of Fig. 3 where the ring is effectively disconnected. This allows us to identify the ratio Λ 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 Fig.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 Tonks-Girardeau 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.

IV. MOMENTUM DISTRIBUTIONS
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 in more detail the actual qubit state, which consists of a macroscopic superposition of countercirculating current states, and assess its visibility through the study of TOF images.
The momentum distribution, which is experimentally accessible in cold atoms experiments via TOF measurements [51], allows to detect the state of current circulation along the ring [52][53][54]. It is defined as the Fourier transform with respect to the relative coordinate of the one-body density matrix ρ (1) where x and x ′ denote the position of two points along the ring's circumference. Although, in general, 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 k's. To adapt Eq. (3) to our lattice system, we is the Wannier function localized on the j-th lattice site, and x j denotes the position of the j-th lattice site. Thereby, Eq. (3) can be recast into wherew(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 Sec. III B for a more detailed discussion).
In the absence of any barrier, the momentum distribution exhibits a perfect ring shape reflecting the conservation of angular momentum, as discussed in Appendix A. As soon as Λ ≠ 0 and close to the frustration point Ω = π, the system enters a superposition of states with different angular momenta, which becomes visible in the TOF images (see, e.g., Fig. 7).
To analyze these images in more detail, 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 Eq. (B4) in Appendix B for the derivation) where the Bessel functions J n correspond to states with angular momentum n, and γ k is the angle along the ring; the parameter ϕ is a function of the flux and the barrier strength (see Eq. (B3) 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 0 and J 1 interfere, giving rise to a peak at zero k and a fringe with ring-shaped symmetry. The visibility this feature increases with the barrier strength Λ. 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 Fig. 7 that, for sufficiently weak interactions, an angular modulation of the ring-shaped momen- tum 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 k values.
The very different TOF images in the regimes of weak and strong interactions can be well understood by recalling the different nature of the superposition state in the various interaction regimes: At zero or weak interactions, the superposition is well described by the so-called NOON state N, 0⟩ + 0, N ⟩, i.e., a macroscopic superposition of states where, e.g., all bosons occupy either the state with zero angular momentum or the one carrying one quantum of angular momentum. The fidelity of this state worsens at increasing interactions [26]. At strong interactions, the system wavefunction is closer to a macroscopic superposition of Fermi spheres [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 Λ c for disconnecting the ring, as identified in Fig. 3. Globally, we observe that good-quality TOF images allowing to easily identify the superposition of current states as a mod- ulated ring structure are found for a ratio Λ U in the vicinity or above the critical line of Fig. 3 (e).
To quantify the visibility of the superposition state in the momentum distribution for different barrier and interaction strengths, we define reflecting the modification in the integrated momentum distribution due to superposition of states induced by the barrier. We find that η is non-monotonic upon increasing the interactions between the particles, while keeping fixed the barrier strength- Fig. 8. This is an effect of the non-monotonic screening of the barrier as a function of interaction strength, first predicted in Ref. [44] through the study of the persistentcurrent amplitude. Finally, we comment on the expected behavior for a system with filling larger than one. In Fig. 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 Eq. (C3). In both cases we note that the TOF images are qualitatively the same as the ones shown in Fig. 7, and therefore our analysis is relevant also for systems with larger number of particles, (which is the realistic case of the experimental realizations). We note that at higher filling, a larger barrier strength is needed with respect to the lower filling case, in order to produce the same superposition and to observe the same TOF.

V. CONCLUSIONS
We considered a system of bosonic atoms loaded in a ringshaped 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 [24]. 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 fully characterize a qubit made of a macroscopic 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 [24], see Figs. 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 Ref. [24] for details).
We studied the resolution and visibility of the qubit, by following two routes: the systematic 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 ∆E 1 and ∆E 2 with system size. Our results indicates that ∆E 1 is appreciable for small and mesoscopic systems while it is suppressed in thermodynamic limit (Fig. 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 [25,26], 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 Fig. 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 visible through time-of-flight 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 visibility of the superposition depends on the relative size between barrier and interaction strengths: the barrier makes the visibility increasing, while the interaction strength, screening the barrier, leads to the opposite phenomenon, yielding a non-monotonous behavior of the visibility. By increasing the filling parameter for fixed U and Λ, 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 visibility 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, sustaining macroscopic quantum effects. The qubit dynamics is visible in a wide range of system parameters. We have identified the ratio U Λ 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 Fig. 3e) and its visibility in TOF images (the visibility is found to be best defined around or above the same critical line). This allows us to conclude that the ratio U Λ on the critical line yields an optimal parameter choice, corresponding to the best trade-off between the momentum distribution visibility and the gap resolution.
In the noninteracting 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 = ̵ h 2 2mL 2 , with m being the particle mass, and L the system size, reads where θ ∈ [0, 2π]. The wavefunction for a state with defined angular momentum is a plane wave ψ(θ) = (1 √ 2π)e inθ , where n ∈ Z to satisfy periodic boundary conditions, and the corresponding spectrum is E n = (n−Ω 2π) 2 . The momentum distribution then reads where R = L 2π is the ring radius, we have defined γ as k x = k sin γ, k y = k cos γ, and J n is the n-th order Bessel function of the first kind. For n = 0 the momentum distribution is peaked at k = 0, while for n > 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 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 the corresponding eigenvalues and eigenvectors reads: where δE = E n+1 − E n , and where cos 2 (ϕ 2) = We then write the wavefunction as where ϕ depends on λ and Ω.
The momentum distribution in this case becomes where an interference term, proportional to cos γ, appears between the two states with defined angular momentum, giving rise to a 2π-periodic angular modulation of the ring shape found previously. This behavior is the same found in Fig. 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.
Appendix C: Methods

Infinite interaction limit: Tonks-Girardeau gas
In the limiting case of infinite repulsive contact interaction between the particles (U → ∞), the so-called hardcore bosons or Tonks-Girardeau gas, an exact approach can be pursed to diagonalize Hamiltonian (1). Since multioccupancy of one site is forbidden by the infinite interaction energy, it can be simplified into where the bosonic annihilation and creation operators have the additional on-site constraints b 2

By applying the Jordan-Wigner transformation
where f i (f † i ) are fermionic annihilation (creation) operators, the Hamiltonian (C1) can be mapped into the one for spineless fermions: This Bose-Fermi mapping is the analogous, for a discrete system, of the one introduced by Girardeau for a continuous system [56]. Hamiltonians (C1) and (C2) have the same spectrum, but non-trivial differences appear in the off-diagonal correlation functions: ⟨f † i f j ⟩ vs ⟨b † i b j ⟩, 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 Fig. 3. This 1D peculiar strongly correlated TG phase has been demonstrated in several experiments on bosonic wires [58,59].

Gross-Pitaevskii equation
In the limit of weak interactions, we adopt a mean-field approximation to simplify the many-body Schrödinger equation. This is the Gross-Pitaevskii (GP) equation for the bosons subjected to a lattice potential, in the presence of a gauge field: where Ψ is the condensate wavefunction, µ is the chemical potential, V 0 is the optical-lattice depth, U b is the strength of the localized barrier, modeled as δ(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 threedimensional scattering length a as g 1D = 2 ̵ h 2 a ma 2 ⊥ [60]. The continuous-model barrier strength U b is connected to the discrete-model one Λ by U b = ΛL M . In absence of the lattice potential, an analytical soliton solution for Eq. (C3) has been found in [44]. 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 time-of-flight images (see Appendix B). In the presence of the lattice potential, we solve Eq. (C3) 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.

Exact diagonalization schemes
a. Working in the full Hilbert space The exact diagonalization (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 exact diagonalization, 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]: In Sec. IV, we considered values of the filling: 5 11, 15 11, 24 9. Correspondingly, the Hilbert space dimensions for that fillings are 3003; 3268760; 10518300.
The non-diagonal part of the Bose-Hubbard 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].
b. Working in the truncated Hilbert space To study the system for larger sizes and larger fillings, the exact diagonalization 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 [64]. 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 Then the truncated Hilbert space can be generated in the following 3 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 [65]: 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 Sec. IV, D 6 = 2345553 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 We assume that: if ⟨n i ⟩+5σ i < K, ∀i then error (in calculating expectation values) is smaller than 0.0006%, if ⟨n i ⟩+4σ i < K then error is smaller than 0.006%, if ⟨n i ⟩ + 3σ i < K then error is smaller than 3% and if ⟨n i ⟩ + 2σ i < 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 Fig. 9 is calculated with the scheme detailed above.

DMRG method
The modified Bose-Hubbard model of Eq. (1) can be quite naturally numerically treated by a Density Matrix Renormalization Group (DMRG) approach, i.e., by optimising a Matrix Product State (MPS) representation of the many-body wavefunction [66,67]. A first requirement, as for almost any numerical treatment of bosons, is to truncate the local Hilbert space down to few states, d = n max + 1, with n max the maximum allowed particle occupation per site. The MPS ansatz is, at first sight, well suited for periodic boundary conditions [68] but a practical implementation of an algorithm over it has to face a number of subtleties and numerical instabilities [67], especially in the case of a non-homogeneous system, like here. At a difference to another recent work of ours [44], then, we decided to opt here not for an explicitly periodic MPS but rather to employ a more standard open boundary (OBC) scheme with a trick. Namely, we represented a ring of M sites as two adjacent stripes of length M 2 linked by only two extremal rungs on first and last site.
On one hand, such an approach implies an effective local dimension d eff = d 2 , as well as a larger bond dimension m (i.e., the matrix dimension), and therefore a priori an extra cost in the tensor contractions involved in the optimization process. On the other hand, though, it allows us to rely on a unitary gauge from the left and the right border, with deep computational advantages: i) The optimization problem can be casted into a standard eigenvalue problem, by means of exact contractions scaling as O(m 3 ). This has to be confronted with the situation for explicit PBC's: there one has to face the instabilities of a generalised eigenvalue solver, whose defining operators are moreover obtained exactly with O(m 5 ) operations [68]; approximate, slightly better scaling, strategies for PBC's are also available, but they work best on very long and uniform chains, which is not the case here [69,70].
ii) The preservation of quantum numbers related to Abelian symmetries, such as the U (1) particle number, is simple to implement, boost the computational efficiency [66], and eliminates an uncertainty source by not invoking any chemical potential [44].
iii) The splitting of a two-site (four-legged) optimised tensor into two single-site (three-legged) ones assumes the clear meaning of an entanglement renormalisation (from m d eff to m states), thus giving a quantitative indication of the performed approximation and permitting a dynamical allocation of symmetry sectors inside the tensors [66]. Such features compensate well the extra d eff cost involved in the contractions with respect to single-site optimization.
In our simulations, [71], we chose the local Hilbert dimension up to d eff = 25 (i.e., n max = 4) in the softer core case U = 1: this can be checked "ex post" to be appropriate, by looking at the decay of site occupation probability and confronting it to other approximations incurring in the algorithm. The other main source of numerical uncertainty is given, of course, by the number of states kept in the RG procedure (i.e., the bond dimension of the MPS ansatz): for moderate ring sizes up to 80-100 sites, as considered here, we have seen that m ≃ 200 already provides reliable results.