Fragmented condensation in Bose-Hubbard trimers with tunable tunnelling

We consider a Bose-Hubbard trimer, i.e. an ultracold Bose gas populating three quantum states. The latter can be either different sites of a triple-well potential or three internal states of the atoms. The bosons can tunnel between different states with variable tunnelling strength between two of them. This will allow us to study; i) different geometrical configurations, i.e. from a closed triangle to three aligned wells and ii) a triangular configuration with a $\pi$-phase, i.e. by setting one of the tunnellings negative. By solving the corresponding three-site Bose-Hubbard Hamiltonian we obtain the ground state of the system as a function of the trap topology. We characterise the different ground states by means of the coherence and entanglement properties. For small repulsive interactions, fragmented condensates are found for the $\pi$-phase case. These are found to be robust against small variations of the tunnelling in the small interaction regime. A low-energy effective many-body Hamiltonian restricted to the degenerate manifold provides a compelling description of the $\pi$-phase degeneration and explains the low-energy spectrum as excitations of discrete semifluxon states.


Introduction
It is well known that bosons at sufficiently low temperatures tend to form Bose-Einstein condensates (BECs), which essentially consist on the macroscopic population of a single-particle state [1]. In absence of interactions, the macroscopically occupied state is the lowest energy state of the single-particle Hamiltonian. When the interatomic interactions are taken into account, and for sufficiently large number of atoms, the main effect is a broadening of the single-particle state, which can be accounted for by a mean-field Gross-Pitaevskii description. In the Onsager-Penrose picture, in a BEC there is only one eigenvalue of the one-body density matrix which is of the order of the total number of particles. This is the largest eigenvalue and is termed the "condensed fraction".
In contrast, an interesting scenario appears when, in absence of interactions, the lowest single-particle states are degenerate. In this case, the many-body ground state may get fragmented [2], as naively the atoms have no reason to condense in only one of the degenerate single-particle states. This implies that a finite number of eigenvalues of the single-particle density matrix are of order of the total number of atoms.
In Ref. [2] the authors describe three notable physical examples which produce fragmented condensates, e.g. the highly correlated regime in ultracold gases subjected to synthetic gauge fields, the two-site Hubbard model and the ground state of a spinor condensate in absence of quadractic Zeeman terms. In the first two cases, interactions need to dominate over tunnelling terms in order to get fragmentation. For instance, for atoms in the double-well it is in the strongly repulsive regime that the condensate fragments in two parts, with half the atoms populating each well. By increasing the potential barrier between the wells, the system enters into the Fock regime (interaction energy dominating over the tunnelling) and fragments in two BECs of equal number of particles. This has been realized experimentally [3,4,5,6]. Notably, the quantum many-body correlations present in fragmented states can find applications in the field of quantum metrology to improve precision measurements [5], and hold promise of being useful in near future technological applications [7].
It is desirable to pin down quantum many-body systems which feature fragmentation even at the single-particle level without explicit spatial separation. In this article we consider a minimal system which fulfills this, and which therefore has fragmented ground states both in absence of atom-atom interactions or for small ones. We consider N identical bosons populating three single-particle states. The bosons are assumed to be able to tunnel between the different states. There are two options that are available with current techniques. The first one would be to trap the atoms in a triple-well potential, with fully connected sites as in [8,9], or aligned [10]. In this case, the three quantum states are the three eigenstates of the single-particle Hamiltonian, which are thus spatially localised. A second option is to consider three different internal states of the boson as single-particle states, with the whole cloud being trapped on the same harmonic potential. Transitions among the three internal states can be induced by Figure 1. Schematic depiction of the system under study. We consider a gas of bosons which can populate three different modes, depicted as green balls. Bosons are allowed to tunnel between the modes. The tunnelling between modes 1 and 3 is taken to be tunable through the parameter γ. For γ = 0 the three sites are aligned and the geometry is essentially one dimensional. For γ > 0 we have a triangular configuration, which for γ = 1 is equilateral. In the limit of γ >> 1 the system is similar to a two-mode one. For γ < 0 we have π-phase tunnelling between sites 1 and 3. means of Rabi coupling. In this case, the three quantum states are not spatially localised. These two cases may be referred to as external or internal, three-mode systems. The properties of such triple-well potentials have been studied previously exploring the effect of dipolar interaction in the system [11,12,13,14,15], the many-body properties of the system [10,16] or the melting of vortex states [9].
In the previous studies, the coupling between the different modes is provided by the quantum tunnelling between the spatially localised modes. This has hindered the exploration of the regime we discuss in the present article. Namely, we will consider triangular setups in which one of the tunnelling terms can be taken negative. This means that a particle tunnelling from one site to the next one acquires a phase of π. Thus, we consider cases when a particle acquires either 0 or π phase when tunnelling around the triangle in absence of interactions. The key idea is to consider systems where tunnelling can be detuned [17] and particularly profit from the recent advances in producing phase dependent tunnelling terms. For external modal configurations, an external shaking of the system along one direction effectively results in a dressed tunnelling term whose sign can be switched from positive (standard) values to negative (π-phase tunnelling) [18]. More recently, a deep laser dip in the centre of a junction has been proposed in Ref. [19] to engineer π-phase tunnelling. In the internal case a phase-dependent tunnelling can be obtained as in Ref. [20]. This paper is organised as follows. In Sec. 2, we present the three-mode Bose-Hubbard Hamiltonian, and briefly recall the many-body magnitudes that are used to characterise the entanglement and correlation properties in the system. In Sec. 3, we first discuss the non-interacting case, whose properties are crucial to understand the onset of fragmentation at the many-body level. Sec. 4 is devoted to analyse the role played by interactions. In Sec. 5 we provide a model developed to understand the appeareance of fragmentation in the symmetric π phase. Finally, in Sec. 6 we summarise our conclusions and provide possible implications for future experiments.

Three-mode Bose-Hubbard Hamiltonian
We assume N ultracold bosons populating three quantum states. As discussed above, they can be different sites, e.g. an ultracold gas confined in a triple-well potential, or three internal states of the atoms. For the time being, we will restrict our system to the former case. We consider tunnelling terms between the three sites, and a tunable rate between two of them. This tunable rate will allow us to explore colinear configurations, closed ones, and also configurations with π phase. Besides the tunnelling, the atoms interact via s-wave contact interactions. This interaction is assumed to be the same between all atoms, which is straightforward in the external case, as there is only one kind of atoms, but would require some tuning in the internal case depending on the species. The three-mode Bose-Hubbard (BH) Hamiltonian we consider is thus, whereâ i (â † i ) are the bosonic annihilation (creation) operators for site i fulfilling canonical commutation relations, andn i =â † iâ i is the particle number operator on the i-th site. J is the tunnelling coefficient between sites 1-2, and 2-3, and U is the atom-atom on-site interaction that can be repulsive U > 0 or attractive U < 0. In our study, sites 1 and 3 are always equivalent with respect to site 2, and the tunnelling between sites 1-3 depends on the particular configuration through the parameter γ. Figure 1 shows schematically the triple-well potentials we have addressed. When γ = 0, no tunnelling exists between sites 1 and 3, which corresponds to an aligned triplewell configuration. When 0 < γ < 1, sites 1 and 3 are connected but the tunnelling between 1-2 (and 2-3 as well) is larger. In contrast, when γ = 1, the three sites are fully equivalent with the same tunnelling rate among them, which can be geometrically interpreted as arranged in an equilateral triangular potential. We will go beyond this symmetric configuration, when γ > 1, by increasing the tunnelling rate between sites 1 and 3 with respect to 1-2 (and 2-3), up to γ = 2 and we will consider negative values of γ besides. The last range can be engineered by lattice shaking along the direction of sites 1 and 3 [18].
The Hamiltonian of Eq. (1) can also be reproduced with a three-component (spinorial, isotopic or atomic) BEC mixture trapped on a single harmonic oscillator with suitable Rabi coupling between the levels (which corresponds to the tunnelling terms).
Our study will concentrate mostly on the repulsive interaction case, which is more prone to be experimentally explored with current setups. We will also discuss briefly the attractive interaction case. In the latter even small asymmetries in the external trapping potentials will eventually have a large impact on the properties of the ground state [21]. Thus, we will introduce small symmetry breaking terms, compared to both the tunnelling and interaction, to consider situations closer to experimental ones and also to control numerically the degeneracies in the problem. We have added the biases between sites 1-3 and 2-3, 13 and 23 , respectively. The considered Hamiltonian now reads:Ĥ =Ĥ + 13 (n 1 −n 3 ) + 23 (n 2 −n 3 ) . (2)

Many-body basis
The system is studied through direct diagonalisation of the many-body Hamiltonian for a fixed number of atoms, N . A suitable many-body basis is the Fock one, which labels the number of atoms in each mode, where |vac stands for the vacuum, and N = n 1 + n 2 + n 3 . The elements of the Fock basis can be expressed as a product state |n 1 , n 2 , n 3 = |n 1 ⊗ |n 2 ⊗ |n 3 . A general many-body wavefunction is thus written as |Ψ = N n 1 ,n 2 C n 1 ,n 2 |n 1 , n 2 , n 3 , where C n 1 ,n 2 is the corresponding amplitude of the Fock state |n 1 , n 2 , n 3 , and n 3 = N − (n 1 + n 2 ).

Coherent states
Coherent states are the closest analogs to classical solutions, in the same way wavepackets are the closest quantum analog to classical trajectories. A general coherent state can be constructed by assuming that all N atoms populate the same single-particle The coherent state reads, Since c i ∈ C, this many-body state has 6 parameters to be determined. Properly normalising the single-particle wavefunction and also realising that there is always an arbitrary global phase, the number of free parameters can be reduced to 4. It is worth stressing that a coherent state as the one defined above corresponds to a fully condensed atomic cloud.

Condensed fractions
The fragmentation properties of the ultracold atomic gas [2] can be investigated by means of the eigenvalues of the one-body density matrix. In our system with N bosons and three different modes, the one-body density matrix of a many-body state |Ψ is a 3 × 3 matrix whose elements are, with i, j = 1, 2, 3. Since |Ψ is normalised, it follows that Trρ (1) = 1.
It is interesting to calculate its eigenvectors (natural orbitals), |ψ i , and eigenvalues, p i , with p 1 ≥ p 2 ≥ p 3 ≥ 0, that satisfy p 1 + p 2 + p 3 = 1. Each eigenvalue of the onebody density matrix gives the relative occupation number of the corresponding natural orbital: p i = N i /N . In a singly condensed system, there is only one large eigenvalue that corresponds to the condensed fraction of the single-particle state, |ψ 1 : p 1 ∼ 1 and N 1 ∼ O(N ), with all the other eigenvalues p j (j = 1) being small ∼ O(1/N ). Instead a fragmented system has more than one large eigenvalue, N i ∼ O(N ), with i = 1, . . . , s, and the rest of eigenvalues p j (j > s) are small ∼ O(1/N ). In this situation the system is not fully condensed, but fragmented, quantum correlations become important and the mean-field approximation fails to describe the system.

Entanglement properties: entanglement entropy and Schmidt gap
Correlations between different subsystems of a many-body quantum system can be quantified performing different bipartite splittings. That is, considering the system as made of two subsystems, tracing out one of the parts, and studying the von Neumann entropy and the entanglement spectrum [22] of the resulting subsystem.
In our case, we consider different spatial partitions of the three-well configuration, e.g. (1,23), (2,13), as in Ref. [14]. From the density matrix of the full system,ρ, correlations between mode i and the rest can be determined by first taking the partial trace ofρ over the Fock-state basis of the other modes. This yields the reduced density matrix on subsystem i,ρ i , that describes the state of this subsystem. For instance, in our system, by tracing out sites 2 and 3, a bipartite splitting of the three-mode system is obtained (1,23), and the reduced density matrix on site 1,ρ 1 = Tr 23ρ , is found to be diagonal in the single mode space of N particles (see Eq.
where |k are states of k particles in mode 1. Note, that the reduced density matrix for state 1 is in general a mixture without a well-defined number of particles. The set of eigenvalues {λ 1 k } is called entanglement, or Schmidt spectrum [22,23,24], and the eigenvalues are the Schmidt coefficients. The Schmidt coefficient λ 1 k is in this case directly the probability of finding k particles in site 1 without measuring the number of atoms in sites 2 and 3. The Schmidt spectrum fulfills i λ 1 i = 1, and contains information about the correlations and the entanglement properties of the state in subsystem 1 with respect to the rest of the system. It is worth recalling that a manybody state is entangled when it cannot be written as a product state. In the case of spatially separated modes, the entanglement we are discussing is spatial.
A measure of the entanglement between the two subsystems is already provided by the single-site von Neumann entropy, which can be computed as S 1 = −Tr(ρ 1 logρ 1 ). Noting that in our case the density matrix,ρ 1 , is already diagonal, see Eq. (8), the entropy can be evaluated from the Schmidt coefficients as S 1 = − i λ 1 i log λ 1 i . In the three-mode system, a signature of the entanglement on the Schmidt spectrum is the following one: if site 1 is not entangled with sites 2 and 3, it should be pure after tracing those sites out, and then the entanglement spectrum would be λ 1 1 = 1 and λ 1 i = 0 for i = 2, 3, . . .. This actually implies a zero of the corresponding von Neumann entropy. A remarkable magnitude defined from the set of λ's is the so-called Schmidt gap, defined as the difference between the two largest and more relevant Schmidt coefficients in the entanglement spectrum of subsystem i, In the case of no entanglement between the subsystems, the Schmidt gap takes its maximum value of 1. A vanishing of the Schmidt gap marks large entanglement between the subsystems. As has been recently pointed out in Ref. [14] for dipolar bosons in triple-well potentials, the Schmidt gap is a good tool to distinguish between phase transitions and crossovers.
Note that no relation exists between fragmentation and entanglement, in the sense that a system can be entangled and fragmented independently. For instance, a system where all the bosons occupy the same spatial mode is nor fragmented neither spatially entangled. A bosonic Josephson junction in the strong repulsive interaction (Fock) regime is a clear example of a non entangled but fragmented state: |N/2, N/2 . In contrast, in the non-interacting regime the bosonic Josephson junction is in a fully condensed state, which has large entanglement between the two sites as seen by the Schmidt coefficients which are λ k = 2 −N N k . Finally, cat states, |Ψ = (|N, 0 + |0, N )/ √ 2, are a well known example of fragmented and entangled many-body systems.

Non-interacting case
In this section we solve the non-interacting case for any value of the tunable tunnelling link γJ. In our calculation we will keep the tunnelling parameter fixed J/h = 1 Hz, since it essentially sets the overall energy scale.
The triple-well configuration we have chosen is intended to remark first, the role played by the topology of the configuration, going from a disconnected chain (γ = 0) to ( 1 ,0 ,-1 ) Single-particle spectrum of the triple-well system as a function of γ. In black, red and green solid lines we depict the eigenvalues E 1 , E 2 and E 3 , from Eqs. (10). The corresponding eigenstates, Eq. (11), are written explicitly for several relevant values of γ. The dashed lines correspond to the spectrum computed with a bias term an essentially double-well system at γ → ∞, through a connected equilateral triangle at γ = 1. It is worth emphasising that we do so by varying the tunnelling strength between sites 1 and 3 and not by, for instance, adding sizeable bias terms to the Hamiltonian, which would indeed also break the symmetry between the three sites. The second important point, is that we consider π-phase tunnelling between the two wing sites (γ < 0), which will indeed have dramatic consequences on the many-body properties of the system.

Explicit symmetry breaking, effect of the bias
The explicit symmetry breaking induced by bias terms in the Hamiltonian (2) breaks the degeneracies present in the single-particle spectrum. As seen in Fig. 2, the crossings at γ = −1 and γ = 1, which occur in the non-interacting system for the ground state and excited states, respectively, become now avoided crossings. For the case we are interested in, when the bias is much smaller than the tunnelling, we can obtain the states dressed by the bias at the degeneracy points by means of first order perturbation theory. In the case of 13 > 0 and 23 = 0, the ground-state manifold at γ = −1 splits, ∆E = 13 , and the corresponding dressed states are, The first one is the new ground state of the system, which has a slightly larger population of particles in site 3. This result is reasonable, since the bias we have considered, 13 > 0, promotes site 3.

Quantum many-body properties of the system
We consider now the effect of repulsive interactions between atoms. We calculate the ground state, by exact diagonalisation of the Hamiltonian (1), for different values of the tunnelling rate γ. In our numerics we will consider up to N = 48 particles. The ground state of the system is characterised by means of, a) coherence properties and fragmentation, b) analysis of the energy spectrum, and c) entanglement spectrum and entanglement entropy.

Analytic results in the |γ| 1 limit
The structure of the single-particle spectrum, see Fig. 2, allows one to build a simple model for |γ| 1. This simple model will serve as guidance to understand many of the properties of the many-body ground state which will be discussed in the forthcoming sections.
We can distinguish two distinct regimes, a) γ 1 and b) −γ 1. In both cases, for sufficiently large values of |γ| the single-particle spectrum approaches E 0 = −|γ|J, E 1 = 0, and E 2 = |γ|J. In the (a) case, the non-interacting single-particle states are |ψ 0 = 1/ √ 2(|1, 0, 0 + |0, 0, 1 ), |ψ 1 = |0, 1, 0 and |ψ 2 = 1/ √ 2(|1, 0, 0 − |0, 0, 1 ), respectively. In (b) the single-particle states exchange their role, with |ψ 2 being the ground state and |ψ 0 the more energetic. For interaction strengths such that N U/J |γ|, the effect of the interaction can be considered perturbatively. In this limit we can ignore particle-hole excitations to the highest single-particle state, and write configurations with k excited atoms as, where N is a normalisation constant andâ † ψ i creates a particle in the single-particle state |ψ i . The energy for this state, up to constant terms, reads, In the non-interacting case, the energy is minimal for k = 0, as expected. For non-zero interactions, this expression has a minimum for where the Int[. . .] stands for the integer part of the expression. Incidentally, for large interactions, N U/J 1 (with |γ|J U ), the minimal energy is obtained for k = N/3, as expected for a Mott insulator of filling N/3. As we decrease interactions, the system goes step by step to k = 0. In particular, many-body states with k and k + 1 excited atoms degenerate if, In this limit, it is quite reasonable to expect that as we increase the atom-atom interactions keeping γ fixed, the system tends to minimise the number of pairs per site, which is achieved by equipopulating the three sites. Eq. (18) predicts ground state degeneracies, i.e. energy crossings, for certain values of |γ| and N U/J. This model works reasonably both for γ −1 and γ 1, where already the structure of the single-particle spectrum starts to resemble the asymptotic one. There is one important difference between γ −1 and γ 1. In the former case, the lowest energy single-particle state is independent of γ (see Fig. 2) and has no population of site 2. This makes the model outlined above fairly accurate to describe the different transitions in the many-body ground state. In the latter however, the lowest singleparticle state only asymptotically approaches the state (1/ √ 2)(|1, 0, 0 + |0, 0, 1 ). In this case, the model only provides a qualitative picture at smaller values of |γ|. Note also, that for |γ| 1, the value of k which minimises the energy, corresponds in this limit to the average population of site 2 in the ground state, while the average population of sites 1 and 3 is (N − k)/2.

Coherence and one-body density matrix
To have an overall picture of the different regimes that we will encounter for different values of N U/J and γ we first study the average populations of the three sites, Ψ GS |n i |Ψ GS /N . Notice that in our system, without bias terms, sites 1 and 3 are equivalent. In Fig. 3 we present results for −2 ≤ γ ≤ 2 and relatively small values of the interaction N U/J ≤ 10, which do not reach the Fock regime. For very small values of the dimensionless parameter N U/J, one does not expect substantial changes from the single-particle case. Indeed in the ground states for γ > −1 the three modes are substantially populated, whereas for γ < −1 the second mode is clearly less occupied than the other two as reflected in the dark region on the left corner of Fig. 3 (b).
In the non-interacting limit, U = 0, the problem becomes a single-particle one and in the case of a symmetric configuration, γ = 1, the many-body ground state can be written as [9,25], in which the average population on each site is, for symmetry reasons, N/3. For large repulsive atom-atom interactions, regardless of the value of γ, the system will fragment in an effort to diminish the number of pairs inside each site, in analogy to, for instance, the double-well [2]. In the large interaction limit, that is in the Fock regime regime: |N/3 ± 1, N/3, N/3 , |N/3, N/3 ± 1, N/3 , and |N/3, N/3, N/3 ± 1 , where the plus (minus) sign refers to a single particle (hole) delocalization.
In the Josephson regime, defined as N U/J 1, the ground state of the system is mostly condensed for γ −0.5. This is reflected in the eigenvalues of the one-body density matrix, see Fig. 4 and left panel in Fig. 5. As already pointed out in the singleparticle spectrum, a very interesting feature is readily found in the vicinity of γ = −1. In this case, the ground state of the system is fragmented in two pieces even in the non-interacting case. As the interaction increases, but still in the Josephson regime, the system is seen to remain bifragmented. From Fig. 5, one can see that the region in the γ-space, around γ = −1, which corresponds to a fragmented condensate, slightly increases with the interaction.
As discussed above, in the limit of the strong Fock regime, the ground state should be essentially fragmented in three pieces (in our case N = 48, which is multiple of 3), corresponding to the ground state |N/3, N/3, N/3 . This is clearly seen in Fig. 6, which depicts the behavior of the three eigenvalues of the one-body density matrix in the γ = −1 configuration for different values of the interaction. In the non-interacting case, the system is fragmented in two condensates (p 1 = p 2 = 0.5 and p 3 = 0), whereas as N U/J increases, p 1 = p 2 < 0.5 decrease and p 3 > 0 increases, fulfilling p 1 + p 2 + p 3 = 1. Moreover, one can see that p i → 1/3 asymptotically. Note however that the origin of bifragmentation is directly related to the degeneracy at the single-particle level and remains in the presence of tunnelling. Trifragmentation requires strong interactions such that essentially tunnelling plays no role and the system can be regarded as three independent condensates. In the right panel of Fig. 5 the behavior of the two largest eigenvalues of the one-body density matrix is shown, as a function of N U/J, for a fixed γ configuration. At γ = −0.99 there is a sharp transition from a singly condensed (at U = 0) to a bifragmented system already for very small values N U/J > 0. However, as γ departs further from −1 (γ = −0.7 and −1.3) the transition between both regimes becomes smoother, and the system remains fully condensed for a larger range of interactions.

Energy spectrum
The model presented in Sec. 4.1 predicts a number of degeneracies in the many-body energy spectrum of the system. Indeed, Eq. (18) predicts N/3 energy crossings for both γ < −1 and γ > 1. As explained above, these predictions are expected to hold for γ −1 and are indeed found in the exact many-body spectrum as seen in Fig. 7. Varying the value of k in Eq. (18) from 0 to N/3 one obtains the large |γ| behavior of the different lines of zero gap in Fig. 7 (a). Similarly, albeit not shown in the figure, for γ 1 the corresponding gapless lines are also found in the system. For γ > 0, also interesting is the closing of the energy gap between the first and second excited states for γ = 1 (see the vertical line in the right panel of Figs. 7). It corresponds to the degeneracy between vortex and antivortex states studied in Ref. [9], whose wavefunctions have been previously obtained in the non-interacting case, see Eq. (12). As the interaction is increased the degeneracy between the corresponding two states remains.

Entanglement spectrum and entanglement entropy
The many-body ground state has been found to be fragmented in the vicinity of γ = −1 and mostly condensed otherwise, for N U/J 1. Besides fragmentation, the three-site configuration considered allows one to study the onset of entanglement and correlations among the three different sites depending on the specific values of the parameters. To characterise the spatial entanglement properties we will use the Schmidt gap and the entanglement von Neumann entropy defined in Sect. 2.4. Due to the structure of our system, with three sites where two of them, 1 and 3, are essentially equivalent, one can consider two bipartite splittings. The first one corresponds to subsystem 1 where sites 2 and 3 have been traced off (1,23). And second, subsystem 2 where subsystems 1 and 3 have been traced off (2,13).
The Schmidt gap corresponding to both bipartite splittings in the (γ, N U/J) diagram is presented in Fig. 8. The first notable feature is that in both bipartite splittings we observe two fan-like radial structures, with straight lines converging to (−1, 0) and (1, c) (with c a constant to be determined later). This structure, similar to the one discussed in Ref. [14], represents the crossings of the first two values of the entanglement spectrum. To better understand the structure, in Fig. 9 we plot the entanglement spectrum (N + 1 coefficients) for both bipartite splittings for a fixed value of N U/J = 2. For γ < −1, each zero of the Schmidt gap shown in Fig. 8 implies a variation of one unit of the most probable population of the untraced mode. This can be understood from the expression in Eq. (A.2), where the reduced density matrix is shown to be diagonal in the Fock basis of the untraced mode.
The simple large |γ| model explained in section 4.1 also explains qualitatively the observed behaviour. The degeneracies in the energy spectrum obtained there correspond to many-body ground states in which the most likely value of the population in site 2 goes from 0 to N /3 and the corresponding one of sites 1 and 3 from N/2 to N/3. We have considered N = 48, see case to N/3 = 16 at γ = −1. As explained in the previous section, the perturbative model captures the physics also for γ > 1. As we can see in the figure, at γ = 1 the most likely population of site 1 is again N/3 = 16 and it keeps increasing as γ is increased.
These crossings described in the entanglement spectrum, clearly visible in Fig. 9, are one of the main signatures of a crossover between two phases, whereas in a quantum phase transition all the Schmidt coefficients degenerate to the same value in the thermodynamic limit (N → ∞) [23]. In our case, it is clear that at (γ = −1, N U/J = 0) all the Schmidt coefficients degenerate even for finite N . This is due to the fact that we already have a single-particle degeneracy.
Certain limiting cases are easily interpreted. For instance, for γ −1, mode 2 is essentially unpopulated and decoupled from modes 1 and 3. This reflects in a singly populated entanglement spectrum for N U/J = 2 in Fig. 9 (lower panel). In this regime, the system was shown to be condensed in Fig. 5, on the single-particle state (|1, 0, 0 +|0, 0, 1 )/ √ 2, which spatially entangles modes 1 and 3, as seen in Fig. 9 (upper panel). In the γ 1 limit the situation is exactly the same, but as explained above, this limit is achieved in practice for much larger values of |γ| than in the γ < 0 case. In the vicinity of γ = 0 the three single-particle states quasidegenerate. This makes that for relatively low interactions, as in Fig. 9, the many-body ground state starts populating the Fock states around N/3 approaching the Mott insulator phase.
These features, described on the full entanglement spectrum, reflect directly on the corresponding von Neumann entropies, see Fig. 10. For instance, the fact that the system approaches the Mott regime for relatively low values of the interaction in the vicinity of γ = 0 reflects almost zero value of the von Neumann entropy for mode 1, which starts to decouple from the other modes. For γ < −1 and small values of the interaction, the system empties mode 2 and decouples it from the other two modes, see Fig. 10 (right panel).

Quantum phase transition for attractive interactions
Similarly to the two-well system, a quantum phase transition can be expected for γ = 1 for attractive interactions [29,30]. In this case, the three sites are completely equivalent. Increasing the attractive atom-atom interactions, the system will minimise energy by clustering the atoms in a single site. In absence of any spatial bias, the ground state of the system will approach the Schrödinger cat-like state, Analogously to the two-site case, this many-body state is not gapped and is quasidegenerate with its first two excitations. In the thermodynamic limit, the transition between the non-interacting state, Eq. (19), and the cat-like state Eq. (20) goes through a transition point at a finite value of N U/J. This transition reflects in the behaviour of the Schmidt gap of the system. As seen in Fig. 8, several zero-Schmidt-gap straight lines tend to converge on a point in the attractive interaction regime on the γ = 1 line. In Fig. 11 we extend the range of parameters to the attractive region and certainly the lines seem to converge at γ = 1 and −9/2 < N U/J < −4, which is where the authors of Ref. [13] predicted the existence of a quantum phase transition. Notably, in contrast with the former phase transition described at γ = −1, this phase transition is only present in the thermodynamic limit.

Detailed analysis of the γ = −1 case
As already shown in Figures 5 and 6, the γ = −1 ground state is found to be bifragmented, with p 1 = p 2 = 1/2, in the non-interacting limit (N U/J = 0). To a good approximation it remains bifragmented for finite but small (N U/J 1) interactions. As N U/J increases further the system approaches the Fock regime and the ground state tends to the well known trifragmented configuration with p 1 = p 2 = p 3 = 1/3. We will here instead focus on the description of the special bifragmented states obtained when N U/J is non-vanishing but small.
At γ = −1 the single particle spectrum, see theH U Hamiltonian becomes, are the eigenvectors ofH U , with an obvious expression for the eigenvalues. Since N hv1 + N hv2 = N , it will be more clarifying to change notations to N hv1 = k and N hv2 = N − k. Then the spectrum becomes which is degenerate in pairs, (k, N − k), except the topmost energy when N is odd. In particular, the ground state is two-fold degenerate. This approximation turns out to be very accurate for small interactions. The effect of including the |e manifold, breaks the degeneracy but does not promote any of the states. The lowest many body states are then well approximated by, with k = 0, . . . , N where the ± sign labels the two states which would be degenerate in energy in absence of coupling with the |e 's. The Ψ (±) k are obviously bifragmented, and have p 1 = p 2 = 1/2. Furthermore, the ground state is a cat-state with all atoms in one or the other of the modes in Eq. (26). These a † hv can be considered creation operators of discrete semifluxon states [32] because, combining Eqs. (22) and Eq. (26), they can be written aŝ A semifluxon is a quantum state which carries half the quantum flux of the vortex states of Eq. (12). In our case, we have a discrete version, as going around the triangle the quantum phase grows from 0 to π, with the phase jump of π imposed by γ = −1. In Fig. 13 we depict the phase structure of the discrete semifluxon state, compared to the usual vortex, Eq. (12) (referred to as "fluxon").
In the full three-mode space, as said above, the degeneracy of each doublet Ψ (±) k , splits and therefore the ground state should be a NOON state of both semifluxon states. In Fig. 12 we depict the overlap between the exact ground state and first excited state and the corresponding k = 0, ± states in Eq. (30). The model is seen to be fairly accurate up to N U/J 10. In an experimental measurement, one should thus find a fractional flux in one or the other direction.

Summary and conclusions
A simple configuration consisting of three sites with a single tunable link has been shown to exhibit a large variety of quantum many-body properties. An important novelty of the proposed scheme is that we have considered cases in which the tunnelling rate in one of the links is either zero-phase or π-phase tunnelling. In both cases, we have performed diagonalisations for finite number of particles, up to N = 48, and have examined the many-body properties of the systems as a function of the atom-atom interaction and the tunable tunnelling rates.
By varying the tunnelling in one link, by means of the parameter γ, which is within reach experimentally as explained in the introduction, one explores configurations ranging from the colinear one γ = 0, the fully symmetric one γ = 1, to the symmetric π-phase one, γ = −1. In the first two cases, the many-body ground state for small interactions N U/J ∼ 1 is highly condensed, with one of the three eigenvalues of the onebody density matrix clearly scaling with the total number of particles. As we approach the γ = −1 point the system departs from condensation and becomes bifragmented for small interactions. This is partly a consequence of the degeneracy in the single-particle spectrum for γ = −1.
Two quantum phase transitions are present in this system. The first one takes place in the symmetric configuration γ = 1, but only for attractive interactions, which therefore makes it difficult to explore experimentally. This phase transition is reminiscent of the one present in other few-site models, e.g. bosonic Josephson junctions, and marks the transition to Schrödinger cat-like ground states in the spectrum [29,30,31]. Interestingly, for γ = −1 we find a second quantum phase transition, which takes place in absence of interactions and which has clear consequences for small repulsive atom-atom interactions. This phase transition has been characterised by the behaviour of the entanglement spectrum for the two possible independent bipartite splittings as we approach the transition point [23]. At the transition the entanglement spectrum degenerates and the corresponding von Neumann entropies exhibit a maximum. Interactions do not wash out the main features of this transition, which becomes essentially a crossover, which has clear consequences at finite interaction, such as the bifragmentation discussed above.
Finally, for the γ = −1 case and small interactions, we have been able to derive an approximate many-body Hamiltonian, which describes the low-energy excitations as excitations of semifluxons. The bifragmentation is readily explained in this limit, and the ground state of the system is found to be a macroscopic Schrödinger cat state of discrete semifluxons with opposite currents.