Robustness of discrete semifluxons in closed Bose-Hubbard chains

We present the properties of the ground state and low-energy excitations of Bose-Hubbard chains with a geometry that varies from open to closed and with a tunable twisted link. In the vicinity of the symmetric $\pi-$flux case the system behaves as an interacting gas of discrete semifluxons for finite chains and interactions in the Josephson regime. The energy spectrum of the system is studied by direct diagonalization and by solving the corresponding Bogoliubov--de Gennes equations. The atom-atom interactions are found to enhance the presence of strongly correlated macroscopic superpositions of semifluxons.


Introduction
One-dimensional chains have attracted a great deal of interest in the past due to their simple analytical treatment both for bosons [1] and fermions [2]. The interplay between superfluidity and the effect of interactions in a one-dimensional system is particularly involved with some notable phenomena depending strongly on dimensionality, e.g. Tonks-Girardeau physics [3,4]. Superfluid properties have been proved in dragged one-dimensional quantum fluids in open geometries [5], together with its breakdown due to interactions [6]. In closed periodic geometries, instead, condensate dragging can lead to a persistent current, experimentally observed in Refs. [7,8,9,10]. Such persistent currents, which have already been experimentally produced in superconducting devices [11,12,13] as well as in polariton condensates [14], opened promising lines of research for future quantum computers [11]. In particular, states with half a quantum circulation or semifluxon states, have been detected in polariton spinor Bose-Einstein condensates [15], and have also been theoretically proposed in superconducting loops with Josephson junctions [16,17].
The physics of persistent currents [9,18] is intimately linked to phase slip events [19,20], which have been widely discussed in (atomic [10,9], superconductor [21,22] and helium [23]) superfluids. They lead to the appearance of objects with topological defects (vortices, solitons) [24,25] that can provoke counterflow superposition [25,26]. Persistent currents may also be produced in ultracold atoms trapped in a few sites of an optical lattice [26,27,28,29], such that one can profit from the nice properties of ultracold atomic experiments, namely the isolation from the environment and the control over the interactions and geometry [30,31]. In this way we have previously shown how topological defects can be produced by manipulating the phase-dependence of a single link [32]. There, a key ingredient was added to the Bose-Hubbard trimer: a single tunable tunnelling link between two modes. This tunable hopping rate can be eventually turned to negative, which in the symmetric configuration induces a π flux through the trimer, thus producing a two-fold degeneracy in the ground state of the single-particle spectrum. These two-fold degenerate states can be interpreted as discrete semifluxon states and provide the basis for the description of the ground state of the system, which is, under certain conditions, a cat state of semifluxon-antisemifluxon states. This symmetric π−flux case is gauge equivalent to a rotating trimer with a π/3 phase change between all sites, a configuration which has received attention in the past [33]. It is also gauge equivalent to the melting of vortex states studied in a higher part of the spectrum in Ref. [34]. Our symmetric configuration provides a specific gauge, which implies a feasible experimental way of producing superpositions of semifluxon states.
In this article, we generalize our earlier studies of three-site configurations to general closed Bose-Hubbard chains with any number of sites. Forcing one of the links to flip the quantum phase, macroscopic superpositions of discrete semifluxon states are found to form the degenerate set of ground states of the system for small but finite atom-atom interactions. This feature makes this setup particularly appealing in contrast to the fully symmetric chains considered before [27], where persistent currents in closed Bose-Hubbard configurations were studied for excited states. There are nowadays a variety of experimental setups capable to simulate such chains, e.g. ultracold atomic gases trapped in a circular array [35]. Beyond ultracold atomic systems, Bose-Hubbard Hamiltonians have recently been engineered in experiments with coupled non-linear optical resonators [36], where a two-mode Bose-Hubbard dimer has been produced with all relevant parameters externally controlled. The extension of these setups to three or more modes is a promising open line of research. Exciton polaritons provide another experimental root to engineer these Hamiltonians, since soon after the twomodes case [37] discrete ring optical condensates have also been reported [38].
Our setup demands a good control on the tunnelling rate between two sites of a Bose-Hubbard chain, both in strength and phase. Ultracold atomic gases have provided several proposals to achieve such dependence of the tunnelling terms in the case of external modes [39,40,41]. Another option would be to replace the external sites for internal ones, building the connected Bose-Hubbard chain of internal atomic sublevels. In this case, the real one dimensional system is replaced by an extra dimension built from the internal sublevels, in the language of Ref. [42]. In this case a phase-dependent tunnelling can be obtained through the Jaksch-Zoller mechanism [43].
The paper is organized as follows. Section 2 introduces the Hamiltonian that describes the Bose-Hubbard chain with a tunable tunnelling. In Sect. 3, we will analyze in detail the properties of the single-particle problem. Section 4 is devoted to study the coherent ground state of the system and excitations over mean-field states. In Sect. 5 we present the macroscopic superposition of semifluxon states that appears for a given set of parameters, and we investigate their robustness in Sect. 6. We summarize our results and provide future perspectives in Sect. 7.

Tunable Bose-Hubbard Hamiltonian
We consider N ultracold interacting bosons populating M quantum states (sites). Following standard procedures [31] the system is described in the lower-band approximation by the Bose-Hubbard (BH) HamiltonianĤ =T +Û, which, under the notation of Ref. [27] takes the form: whereâ(λ) (â † (λ)) are the bosonic annihilation (creation) operators for site λ, fulfilling canonical commutation relations. J is the tunnelling amplitude between consecutive sites and the parameter U takes into account the on-site atom-atom interaction, which is proportional to the s-wave scattering length and is assumed to be repulsive, U > 0. Attractive atom-atom interactions can also be produced and would more naturally lead to macroscopic superposition states but they are fragile against instabilities. We will analyse the ground state structure and the excitation spectrum as a function of the dimensionless parameter Λ ≡ N U/J, which gives the ratio between interaction and tunnelling rate. As mentioned in the introduction, a crucial ingredient in our model is the presence of a single tunable link. In practice, the tunnelling between sites 1 and M can be varied through a parameter γ, which is taken to be real. This allows to study very different configurations, for instance, an open 1D Bose-Hubbard chain when γ = 0, and a symmetric non-twisted (twisted) closed chain when γ = 1(γ = −1). In the latter two casesĤ is particularly simple. It can be shown that, for the special cases γ = ±1,Ĥ is gauge equivalent to a symmetric Hamiltonian (equal couplings) with a total flux Φ = 0 and Φ = π, respectively: In the previous equation, the site λ = M + 1 corresponds to site 1 (periodicity of the chain). This many-body Hamiltonian describes a Bose-Hubbard 1D circular chain, rotating around its symmetry axis with angular frequency Φ [33,27,44]. Therefore, the physics of the special γ = ±1 cases will be essentially equivalent to the ones considered prior. However, it is worth emphasizing that such local gauge transformation relating the Hamiltonians in Eq. (1) and Eq. (2) does not exist for a general value of γ.

The non-interacting problem
3.1. Single-particle problem For a single particle the problem reduces to finding the eigenvalues µ and eigenvectors χ(λ), with λ = 1, . . . , M , ofT : Using standard techniques in solving tight binding problems, one finds that the solutions can be either even (symmetric): χ (S) (λ) = χ (S) (M + 1 − λ) or odd (antisymmetric): . They can be written conveniently in terms of Bloch phases φ as where the normalization factor of each wave function is N ± = (M ± sin(M φ)/ sin φ)/2. The respective Bloch phases satisfy the implicit equations, for the even solutions, and for the odd ones. In terms of the Bloch phases, the eigenvalues (µ) are then given by For arbitrary γ, all the real solutions that can be extracted from those figures fulfill the condition 0 ≤ φ ≤ π. Therefore the corresponding energies, Eq. (7), are in the band −2J ≤ µ ≤ 2J, and the associated eigenstates are "bulk states", see discussion in Sect. 6. This is in agreement with the results in Fig. 2, which shows the single-particle spectrum as a function of γ for M = 3, 4, 5 and 6. However, there are some special states whose energy is not contained within the previous interval. They correspond to "surface states", and will be discussed later.
The poles of Fig. 1 (which give the solutions when |γ| → ∞) are determined by the zeros of the denominators in Eqs. (5) k,(|γ|→∞) = πk/(M − 1) where k accounts for the odd (even, including k = 0) natural numbers for the symmetric (antisymmetric) wave function. As seen in Fig. 1 each φ corresponding to a finite γ is bounded by those of two consecutive poles φ (S,A) k,(|γ|→∞) , where φ −1,(|γ|→∞) has to be taken as 0. These inequalities show that there are no crossings of single-particle levels of the same parity. This feature can be readily seen in Fig. 2, since even single-particle levels (filled symbols) do not intersect levels corresponding to other even-parity states, independently of the number of sites (and the same occurs for odd-parity states, represented by open symbols). Note also that the curves associated to the solid red (dashed blue) lines in Fig. 1 are monotonically decreasing (increasing). Thus a monotonic variation of γ gives also a monotonic variation of its corresponding energy in Fig. 2.
and at small η, the r.h.s. behaves as 1 + M/2 η 2 and leads therefore to γ > 1. For large η the r.h.s. rises as e η and guarantees that there will be solutions for any γ > 1.
For the antisymmetric solutions one has to introduce also φ = iη and the r.h.s of Eq. (6) becomes and for small η in agreement with the lower bound in Fig. 1. The asymptotic behaviour for large η is now γ − exp(η). Figure 1 also shows that there is a similar problem with the solutions associated to states with highest φ, corresponding to the surface states above the grey region in Fig. 2. In this case one has to set φ = π + iη to find the missing solution. In analogy to the expression of γ as a function of the imaginary Bloch phase, the eigenenergies change as well, µ = −2J cosh η, which at large values of η becomes µ = ∓2Jγ, for the even (−) and odd (+) surface states. This explains the asymptotic linear behaviour of the energies of these states shown in Fig. 2. 3.3. The special γ = ±1 cases These expressions imply equidistributed particle population within all the sites but with a phase variation whose gradient between sites lead to an azimutal velocity. The γ = 1 case is the commonly considered situation in the literature, as it also appears in the usual tight-binding models in condensed matter systems. Concerning the currents, there is one important difference between γ = 1 and γ = −1. In the former case, the ground state of the system corresponds to φ = 0, which is a currentless state. It is always non-degenerate and its eigenenergy is µ gs = −2J, independently of the number of sites M . The first two excited states are degenerate and correspond to φ = 2π/M and −2π/M . They are the discrete version of the usual vortex states (also called fluxons) with a circulation of ±2π [34,27]. In contrast, the ground state of a BH closed chain with γ = −1 is always degenerate. It is spanned by the eigenvectors corresponding to q = 0 (φ = π/M ) and q = M − 1 (φ = −π/M ) in Eq. (11), which are discrete semifluxon/antisemifluxon states (half-vortices with ±π circulation, see Appendix A for further discussion of their properties). The energy gap ∆ ≡ µ ex − µ gs between the degenerate ground states and the first excitation is As explained above, the γ = −1 case can be related by a local gauge transformation to a BH Hamiltonian in which each hopping induces a π/M phase, as in Eq. (2). In that case, considered for instance in Ref. [33] for M = 3, the degeneracy takes place between the fully symmetric ground state for γ = 1, and one of the vortex states.

Many-body coherent states
The eigensolutions ofT described in the previous Section define a new basis for the single particle states, the "mode" basis. The associated creation and annihilation operators will be written asb † q andb q , q = 1, . . . , M , so that with all sums running from 1 to M . From unitarity we have and The coherent states, defined as will be also named mean-field states for N bosons. For the special cases when γ = ±1 we will see later that working in the "flow" basis, Eq. (11), can be more convenient. We define the corresponding creation operators as where the coefficients are given in Eq. (11). For γ = −1, the cases q = 0 and q = M − 1 correspond to the semifluxon,b † sf , and antisemifluxon states,b † asf , already discussed.

The mean-field ground state
The expectation value of the Hamiltonian (1) in a coherent state leads to the Gross-Pitaevskii (GP) energy functional, which will be useful to obtain the excitations over mean-field coherent states. The expectation value ofT can be readily computed by using thatâ † |n = √ n + 1 |n + 1 andâ|n = √ n |n − 1 . Since, Therefore, the matrix elements ofT are, One can analogously compute the expectation value of the interaction part of the Hamiltonian to get: Once the matrix elements of the Hamiltonian are obtained, the corresponding GP meanfield equation can be derived by adding a Lagrange multiplier to conserve the norm µ q λ |χ q (λ)| 2 and differentiating with respect to χ * q (λ). We arrive at the mean-field equations, which in matrix form becomes Eq. (3) for the non-interacting case. Let us point out that for a consistent derivation one has to assume that the χ q (λ) in Eq. (13) are the ones that follow from solving the GP equation including the interaction,Û. In general these χ q (λ) will be different from the ones found without interaction.

Elementary excitations
We will now assume that to a good approximation, even when U = 0, the simplest excited states can be approximated as one atom being promoted from the coherent state |Ψ (N ) q made of all the atoms occupying the χ q to an excited orbital χ p , both p and q orbitals being eigenstates of the non-interacting single particle Hamiltonian,  23) and (25), replacing q and p by 0 and 1. Note also that for −1 < γ < 1 the ground and first excited states change symmetry, see Fig. 1. The cross and star symbols are the more accurate Bogoliubov-de Gennes predictions for γ = 1. In all cases, N = 4 and M = 6.
The expectation value of the Hamiltonian (1) is the sum of the expectation value of the kinetic term and the interactions. The former element is: The kinetic energy cost of the promotion of one particle from the q mode coherent ground state to the p mode is: A similar procedure can be followed to obtain the expectation value ofÛ together with the excitation energy cost due to interactions δU ≡ Ψ Figure 3 shows the calculated lowest excitation energy, δE = δT + δU for a range of values of γ and two particular values of U . In addition, at γ = 1 the energy gap predicted by the Bogoliubov-de Gennes approach computed as in Eqs. (20) and (21) of [27] has been added (cross for U = 0.5 J and star for U = 0.1 J). The figure shows that for small values of U and except in the vicinity of γ = −1 the present approximation can explain most of the effect of the interaction. The singular point γ = −1 will be discussed in the next Section.

Two-orbital approximation for the ground state manifold
As explained above, the non-interacting ground state for the γ = −1 configuration has a two-fold degeneracy between the two semifluxon states. In Ref. [32], it has been shown that for M = 3 and small interactions, the low-energy states of the system can be described by a two-mode model involving only these two single-particle states. For a closed chain with M sites, one can expect that the description of the system as a macroscopic superposition of two countercirculating semifluxons can be generalized. Moreover, the single-particle energy gap of Eq. (12) can be expected to protect the persistent currents created in the ground state manifold in ultracold atomic physics experiments. For M 8 the gap is of the order of the tunnelling J. Thus, for few sites and small interactions (N U J), the physics can be restricted to the degenerate ground state manifold.
One can generalize the procedure followed in Ref. [32], by writing the creation and annihilation operators in the coherent flow basis, and truncating such decomposition to the semifluxon (sf) and the antisemifluxon (asf) states, i.e., and the corresponding expression forâ(λ). One can then rewrite theÛ operator by taking into account that sums like M λ=1 exp(i2λπ/M ) vanish, aŝ whereN x =b † xb x with x = sf or x = asf are the number operators of semifluxon or antisemifluxon states. Since N = N sf + N asf is constant, the last term in Eq. (27) corresponds to a global energy shift that can be neglected. The eigenstates ofÛ are the new "Fock" states  In Fig. 4 we plot the energy spectrum obtained by exact diagonalization of the many-body Hamiltonian, Eq. (1). The band structure of the energy spectrum for small interactions, Λ 1, can be understood by means of the number of atoms and the degeneracy of the flow basis elements (see Appendix A). In the figure we show two sets of parameters (M, N, Λ) corresponding to (4, 10, 0.5) and (5,8,2). The spectrum of the first set shows traces of the degeneracy pattern present in the non-interacting case whereas for the second set the gaps close in the middle and slightly in the upper part of the spectrum. The inset in Fig. 4 shows the comparison between the low-lying exact many-body spectrum with the prediction of Eq. (27). This approximation turns out to be very accurate for small interactions Λ < 2. For the two cases considered in the inset one can see that the model provides a good description for Λ = 0.5, but it starts to deviate for larger values, such as Λ = 2.
For small Λ, the eigenvectors belonging to the low energy manifold can be well approximated by using the Fock basis, |N sf , N asf . The integer index k runs from 0 to N , and the ± sign labels the two states, which would be degenerate in energy if Eq. (27) was exact.
In this approximation, |Ψ ± 0 corresponds to the two-fold degenerate ground state. A similar two-orbital approximation and macroscopic superposition of superfluid flow was considered in Ref. [33] for a three well system in a different gauge.

Bogoliubov-de Gennes spectrum
The spectrum of elementary excitations in the weakly interacting regime can be studied within the Bogoliubov-de Gennes (BdG) framework. We follow the model presented in Ref. [27] for a circular array of Bose-Einstein condensates with the same tunnelling rate between the sites (γ = 1). For the case of γ = −1, the equations are formally the same but with φ q = 2π M (q + 1 2 ) instead of φ q = 2π M q (see the discussion below Eq. (7)). The Bogoliubov excitation spectrum is constructed over mean-field states defined as macroscopically occupied modes in the flow basis, Eq. (15), with q = −1, 0, 1, . . . , M − 2 (or equivalently q = 0, 1, . . . , M − 1). In the previous expression the corresponding coefficients of the creation operatorb † q areχ q (λ; γ = −1) defined in Eq. (11). The periodicity of the system imposes that q is a cyclic index, with period M : for example, q = 0 (−1) and M (M − 1) are equivalent.
The excitation energies relative to a macroscopically occupied state can be obtained with the same procedure developed in Ref. [27]: with The cyclic index k runs from 0 to M − 1, and it is interesting to note that when k = 0 (M ), the Bogoliubov spectrum yields E (±) k=0 = 0, which corresponds to remain in the unperturbed mean-field state (without any excitation). Figure 5 shows the exact many-body spectrum and the BdG excitations (indicated by arrows) relative to the mean-field like states of a system of N = 4 atoms in M sites. We have calculated the BdG spectrum in a chain with M = 3 and M = 5 sites  Tables B1 and B2). Moreover, the solutions must fulfill that the excitations relative to the ground state are positive, whereas those relative to the highly excited macroscopically occupied state in the highest band must be negative.

Robustness of superconducting flows when γ = −1
The macroscopic superpositions of semifluxon states are predicted to appear for low interactions in the special γ = −1 case. Our focus here will be on studying the presence of such macroscopic superpositions of superfluid flows when γ = −1. Two indicators will be used to signal the presence of macroscopic cat-like states as the ones written in Eq. (29) with k = 0. The first one is the overlap between the cat states and the exact solutions resulting from the numerical diagonalization. As discussed earlier, for Λ 1 the ground state is almost doubly degenerate, then, {|Ψ + 0 , |Ψ − 0 } is a suitable basis for the lower manifold, as predicted by the two-orbital model. We therefore define the overlap determinant as, with |Ψ  The second indicator is the fragmentation of the ground state of the system, given by the eigenvalues n i of the single-particle density matrix,ρ  Note that the overlap determinant would be strictly zero if the quasi degeneracy is absent. In this case, even though the ground state may still be well described by |Ψ + 0 , we would get a zero overlap between the two manifolds. This is the reason why the overlap becomes abruptly zero in the vicinity of γ = 1 signalling a level crossing in the many-body spectrum between the first and second excited states, see Figs. 7 and 8. This level crossing is directly related to the crossing found at γ = 1 in the single-particle spectrum, see Fig. 2. Even though the two-fold degeneracy is broken, we find that the actual ground state of the system has a sizeable overlap with the |Ψ + 0 state in a broader region of parameters, as seen in Fig. 7(b). There, for N = M = 5 we find that Λ 1 is a good candidate to find counterflow cat states for a broad range of γ. The overlaps between the ground state and the two quasi-degenerate cat states in the vicinity of γ = −1 depend critically on the number of particles. In the case N = M = 5 for γ −1 we find that the overlap between the ground state and |Ψ + 0 goes from almost one for γ −1 to close to zero for γ −1 for Λ 1. While the situation is the opposite for the overlap with |Ψ − 0 . A much more symmetric situation is found for even number of particles, which does not show this change, as seen in Fig. 7(d,e). This behaviour can be understood from the two-orbital description of Sect. 5.1, details are beyond the scope of the present manuscript.
Further increasing the interactions, Λ > 100, the condensed fraction decreases quickly to the expected value 1/M for the Mott insulating regime, Fig. 6. For integer filling, in the Λ → ∞ limit the ground state is non-degenerate and gapped (Mott insulator of the corresponding filling). Thus, it is expected that already for any finite value of Λ the two-fold degeneracy predicted by the two orbital model would be only approximate. For fractional fillings, in contrast, in the large interaction limit the ground state will be degenerate, with the surplus particles delocalized in the chain. In this case the numerics shows that the two-fold degeneracy predicted by our model is present for all values of Λ.
Besides the condensed fraction and the overlap, a crucial feature of a macroscopic superposition is its robustness, characterized by the presence of an energy gap that protects the superposition to be affected by excitations that could involve larger-energy excited states. In Fig. 9 we present the energy difference between the ground state and the first excited state, panel (a), and the energy difference between the latter and the second excited one, panel (b). There are two important properties exhibited by the peculiar case of γ = −1. First, for zero interactions, the ground state has a very large degeneracy, N + 1, stemming from the combinatorial factor, i.e. N particles populating two single-particle states in N + 1 ways, see Appendix A. For very small interactions, Λ 0.5, this is reflected in the very small gap 0.01 J both between the ground and first, see Fig. 9(a), and between the first and second excited many-body eigenstates, see Fig. 9(b). Secondly, for γ = −1 and Λ 5, a gap starts to open above the ground state. This gap opening increases the robustness of the macroscopic superposition, which indicates that experimental observation of these states should be feasible.
When U = 0 the densities can be easily constructed using the wavefunctions given by Eq. (4). When γ > −1, the ground state is symmetric, see Fig. 1, and the density is given by, with φ determined from Eq. (5). These densities are linear in N , have a maximum at λ m = (M + 1)/2 and are site symmetric with respect to λ m . The simplest example is the case γ = 0: then the solution of Eq. (5) is φ = π/(M + 1) and which shows that at the maximum ρ M = 2N/(M + 1) and at sites 1 and M , ρ = 2N/(M +1) sin 2 (π/(M +1)). Analogous expressions, involving hyperbolic functions, can be easily derived for the surface states.
Similar expressions can be written when γ < −1, now using the odd single-particle solutions. The exact results shown in Fig. 10 for Λ = 0.1 are almost identical to these simple predictions, except when γ is very close to −1 where the breaking of the degeneracy of the odd and even solutions when U = 0 has to be taken into account.
The effect of varying γ and Λ on the density of the cloud is also particularly pronounced. For γ = 1, the population is equal within all the sites, as the system is rotationally symmetric. For γ = −1 the situation is different due to the quasidegeneracy at the ground state level. For small values of the interaction, the ground state is well represented by the cat states (29), which have an equal amount of fluxon and antifluxon components resulting also in a constant density along the chain, see Fig. 10. In the Mott regime, the system is equipopulated again, regardless of the value of γ, see the Λ 10 results in Fig. 10. Already for |γ| 1.5, the density approaches the one built from the surface modes described in the first section, with population peaked on the sites around the tunable link. In Fig. 10 we provide a broader picture for filling factor 2. Again for large interactions, both the central density and the density at the extremes approach the N/M = 2 limit. Interestingly, the region of macroscopic superposition of fluxon-antifluxon states reflects in an almost equipopulation of all sites for all values of Λ. As discussed above, away from γ = −1 and for lower interactions, the cat-like state is less robust, resulting in a higher density in the center (extremes) as γ increases (decreases). Monitoring the density of the chain can therefore be a good indicator of the macroscopic superposition states expected. For instance, the chain could be initially prepared at small but nonzero interactions and γ = 1. Turning the tunable link from γ = 1 to γ = −1, the density in the center will grow, reaching very large values for γ −1. At γ = −1 the chain should again be equipopulated. This transition from having almost zero population in the extremes to equipopulation would signal the regime of macroscopic superposition of fluxon-antifluxon states.
The transition from a condensed to a fragmented state consisting in a macroscopic superposition of counterpropagating flows has been described for finite number of particles and sites using exact diagonalization techniques. It has been shown to take place at small but nonzero interactions, which actually helps to stabilize the superposition. Now we will stretch our numerical diagonalization techniques to explore the thermodynamic limit at small Λ. We have used the ARPACK package to go up to Hilbert spaces of dimension 10 6 . This allows us to explore the behaviour of the macroscopic superpositions of semifluxon states with up to N = M = 13, as shown in Fig. 11. To roughly explore the thermodynamic limit we have performed extrapolations using up to quadratic terms in 1/N , shown in the figure with dotted lines. The transition to the macroscopic superposition phase is clearly seen to be much more robust at Λ = 10 than Λ = 0.1. In the former, for instance we find values of the overlap determinant of about 0.25 for γ = −0.5, −1 and −1.5 in the extrapolation to the thermodynamic limit. The corresponding ground state fragmentation is close to 0.5 in agreement with the predicted bifragmented nature of the superposition state. For smaller Λ the superposition is only found at values γ −1: departures from it lead to a quick loss of the bifragmentation and to small values of the overlap determinant.

Summary and conclusions
In this paper we have studied the role of a tunable link in a Bose-Hubbard chain of arbitrary size. In the case of non-interacting closed systems, we find counterpropagating persistent current states in the upper part of the spectrum with γ = 1 [34], and in the ground state for the case of γ = −1, constituting a cat state of macroscopic flows. In this latter case, we have also analyzed the Bogoliubov excitations over mean-field states when interaction is present, by following a procedure similar to the one presented in Ref. [27].
We have analyzed the robustness of these counterflow cat states by studying the energy spectrum, the density profile, the condensed fraction and the overlap of the non-interacting ground state manifold with the ground state computed by exact diagonalization, as a function of the interaction and γ. We have found that macroscopic superpositions of counterflow in weakly interacting Bose gases in a closed Bose-Hubbard chain with γ = −1 is more robust against fluctuation of the parameters, and is protected by a larger energy gap from excited states. Quadratic extrapolation to the large N limit in the filling factor N/M = 1 regime predicts the existence of such superpositions in the thermodynamic limit for small but non-zero interactions, Λ 5.
The production of macroscopic superpositions of semifluxon states in interacting many-body systems opens a new possibility to obtain persistent currents which may be useful in quantum computation and quantum simulation. An important feature of the results described in this article is that they are applicable to rings of any number of sites and thus can be engineered in a large variety of experimental setups. Thus, closed Bose-Hubbard chains with a single tunable link provide a versatile system in which macroscopic superpositions of flow states can be produced with a tunable twisted link, in which the tunnelling can be varied both in strength and phase. Among the possibilities that exist nowadays, cold atoms [35] and coupled non-linear optical resonators [36] are promising quantum devices to produce countersuperflow states.
The band structure of the non-interacting spectrum and also for small interactions Λ < 1, see Figs. 4 and 5, can be understood by means of the number of sites, the number of atoms, and the degeneracy of the flow basis elements.
For example for M = 3 sites and γ = −1 the single-particle spectrum has a ground state with two degenerate eigenvectors, |fl and |afl , and an excited state |sf + fl . For a system with N atoms, the total number of bands is N + 1 that corresponds to the different possibilities of distributing N atoms in a two level system with occupancies N 1/2 ≡ N sf + N asf and N 3/2 ≡ N sf+fl , with the restriction N = N 1/2 + N 3/2 . Inside each band, the number of quasi-degenerate levels can be calculated by counting the number of different possibilities (N sf , N asf ) to set N 1/2 atoms between two degenerate levels |sf and |asf , where N 1/2 ≡ N sf + N asf . The quasi-degeneration is deg 1/2 = N 1/2 + 1.
In tables A1 and A2 we present the distribution in bands and the quasi-degeneracy inside each band, for the many-body spectrum of N weakly-interacting atoms in M sites. Figure 5 shows the band structure of the exact many-body spectrum for N = 4 atoms in M = 3 and M = 5 sites, for small interactions. The band degeneracy predicted in Table A1 and A2 is in agreement with the quasi-degeneracy inside each band obtained numerically (see Fig. 5).  = 4  b1  b2  b3  b4  b5  b6  b7  b8  b9  b10  b11  b12  b13  b14 b15 4  3  3  2  2  1  2  1  0  1  0  1  0  0 Table A2. Distribution of N = 4 atoms in the single-particle flow states of M = 5 sites with γ = −1: |sf , |asf , |sf + fl , |asf + afl , |sf + 2fl . The many-body states in the flow basis are |N sf , N asf , N sf+fl , N asf+afl , N sf+2fl . The corresponding mean-field like states are marked in boldface and correspond to: |4, 0, 0, 0, 0 , |0, 4, 0, 0, 0 , |0, 0, 4, 0, 0 , |0, 0, 0, 4, 0 , and |0, 0, 0, 0, 4 . This case corresponds to a three level system with occupancies N 1/2 , N 3/2 and N 5/2 with N = N 1/2 + N 3/2 + N 5/2 and 15 bands. The two first ones are degenerate and belong to the lowest band, and the latter is the unique state in the highest band. One can consider a small rotational bias to break the degeneracy, such as the many-body ground state of the system is either the macroscopically occupied semifluxon mode |Ψ , there are three possible values of the index k = −1, 0, 1, and therefore two elementary excitations k = ±1. In the BdG framework [27], they can be understood as the building up of a fluxon (vortex-like excitation) when k = 1 or an antifluxon (antivortex-like excitation) when k = −1.
Let us consider the macroscopically occupied semifluxon mode |Ψ (N ) q=0 = |N, 0, 0 , an excitation with k = −1 will lead the system to the excited state |N − 1, 1, 0 which corresponds to the promotion of one particle from |sf to |sf + fl , whereas a k = 1 excitation will lead to the final state |N − 1, 0, 1 . In Fig. 5 the BdG excitations are shown for a system with M = 3, N = 4 and U/J = 0.1. When the macroscopically occupied mode is the semifluxon mode |Ψ (N ) q=0 = |4, 0, 0 , a k = −1 excitation leads the system to the excited state |3, 1, 0 which is the third one in the lowest band, see Table A1 in Appendix A. The BdG excitation energy is E q=0 k=−1 /J = 0.0969 which is in good agreement with the excitation energy calculated by exact diagonalization ∆E q=0 k=−1 /J = 0.0992 (see Table B1). A k = 1 excitation relative to |Ψ (N ) q=0 = |4, 0, 0 , leads to the many-body state |3, 0, 1 , whereas a k = −1 excitation leads to |3, 1, 0 . Both final states are degenerate in the fourth band (see Fig. 5), and therefore the BdG prediction gives the same excitation energies E q=0 k=1 /J = E q=0 k=−1 /J = 3.0969 in quantitative agreement with the exact excitation energy ∆E q=0 k=1 /J = 3.0998 (see Table B1). In Table B1 we compare the excitation energies calculated within the BdG formalism and the ones obtained by exact diagonalization of the BH hamiltonian for M = 3 and N = 4 and U/J = 0.1. As can be seen in the table, the agreement for weakly interacting systems is very good. For N atoms in a BH chain with M = 5 sites the flow basis is |N sf , N asf , N sf+fl , N asf+afl , N sf+2fl . There are 5 macroscopically occupied states (q = −1, 0, 1, 2, 3), see Table A2 for N = 4 atoms. Two of them are degenerate in the lowest band (b1): |Ψ 0 q=0 = |N, 0, 0, 0, 0 , and |Ψ q=2 = |0, 0, 0, 0, N . For a mean-field state, i.e. for a fixed q, the BdG framework will provides 5 possible excitation energies labeled by k = −1, 0, 1, 2, 3 (or an equivalent set of values due to the periodicity).
Again, one can consider a small rotational bias such as the many-body ground state of the system is the macroscopically occupied semifluxon state |Ψ (N ) q=0 = |N, 0, 0, 0, 0 . With another bias the ground state of the system can be |Ψ (N ) q=−1 = |0, N, 0, 0, 0 , that corresponds to the macroscopically occupied antisemifluxon state.
We have calculated the BdG excitations related to the macroscopically occupied states for N = 4 atoms and U/J = 0.1. In Table B2 we compare the BdG spectrum with the excitation energies obtained from the BH Hamiltonian. We have discarded all solutions which do not fulfill the BdG normalization [27]. From this comparison, it follows that the BdG framework provides a good description for low-lying excitation energies when the interactions are small.