Revealing quantum effects in bosonic Josephson junctions: a multi-configuration atomic coherent states approach

The mean-field approach to two-site Bose-Hubbard systems is well established and leads to nonlinear classical equations of motion for the population imbalance and the phase difference. It can, e.g., be based on the representation of the solution of the time-dependent Schrodinger equation either by a single Glauber state or by a single atomic (SU(2)) coherent state [S. Wimberger et al., Phys. Rev. A 103, 023326 (2021)]. We demonstrate that quantum effects beyond the mean-field approximation are easily uncovered if, instead, a multi-configuration ansatz with a few time-dependent SU(2) basis functions is used in the variational principle. For the case of plasma oscillations, the use of just two basis states, whose time-dependent parameters are determined variationally, already gives good qualitative agreement of the phase space dynamics with numerically exact quantum solutions. In order to correctly account for more non-trivial effects, like macroscopic quantum self trapping, moderately more basis states are needed. If one is interested in the onset of spontaneous symmetry breaking, however, a multiplicity of two gives a big improvement towards the exact result already. In any case, the number of variational trajectories needed for good agreement with full quantum results is orders of magnitude smaller than in the semiclassical case, which is based on multiple mean-field trajectories.


I. INTRODUCTION
The Bose-Hubbard (BH) model of S interacting (bosonic) atoms in optical lattices is the basis of many state of the art experimental [1][2][3][4] as well as theoretical efforts [3,5,6]. The cold atom Hubbard tool box introduced in [5] puts a focus on strongly interacting many-body dynamics and embraces the fields of quantum optics, quantum computation and solid state physics. The BH model is a paradigm for the rich physical phenomena exhibited in these areas, such as there are quantum phase transitions between the superfluid and the Mott insulator phase [1], self-trapping in bosonic Josephson junctions [7], and quantum chaology [6], to name just a few.
Restricting the amount of lattice sites makes the quantum dynamics of the BH model easily tractable numerically for moderate particle numbers. Recent theoretical work has thus focused on the cases of four (and six) sites [8] with different levels of approximation: exact, semiclassical and classical (mean-field, or truncated Wigner approximation (TWA)).
Also the trimer (ring) case has been studied, due to the facts that it is leading to the melting of discrete vortices via quantum fluctuations [9] and that it is the smallest system that displays a mixed phase-space mean-field dynamics without an external driving term [10,11]. This system has also been dealt with using a group theoretical [12,13] and a semiclassical time-domain approach [14]. With an additional drive (periodic kicks) even the double well case is showing signatures of chaos [15]. Furthermore, the case of two wells without external driving has been extensively studied. The system dynamics has, e. g., been investigated both in a mean-field classical approximation and fully (and perturbatively) quantum mechanically [7,[16][17][18][19][20], as well as also semiclassically, using a phase space picture [21], or employing the Herman-Kluk propagator [14,22]. This same propagator has also been used in a semiclassical time-domain study of the single well problem [23]. Furthermore, the driven single well problem has served as a model in a study of dynamical tunneling [24].
An important lesson from the vast literature is that semiclassical approaches do well in reproducing the full quantum results, while the mean-field and/or truncated Wigner approach have their limitations. TWA does, e.g., not allow for the investigation of revival phenomena, present in the quantum dynamics [20]. In contrast, the macroscopic quantum self trapping effect in bosonic Josephson junctions could already be uncovered using a meanfield approach based on the Gross-Pitaevskii equation [7]. It turns out that mean-field theory predicts the transition to macroscopic quantum self trapping at too large values of the on-site interaction strength, however [25].
In the following, we will focus on the quantum dynamics in the case of two wells, for which the direct experimental observation of tunneling and self trapping has become possible [26].
Theoretically, this case has been reviewed in [27] as well as in [28], where the exact solubility of the eigenvalue problem in terms of the Bethe-Ansatz has been reviewed. Furthermore, a fresh look on finite size (i.e., finite particle number) effects in the mean-field dynamics of those Josephson junction systems has been given by Wimberger et al. [25]. These authors have used a so-called atomic or SU(2) generalized coherent state [29], to uncover meanfield 1/S corrections to the more familiar mean-field results based on standard Glauber coherent states. We will also employ those favorable number conserving SU(2) states here.
We will not use them in a mean-field spirit, however, where just a single state is taken to solve the time-dependent Schrödinger equation (TDSE). In contrast, we will investigate what happens if we allow for non-trivial multiplicity, which for reasons of simplicity we first choose to be just two, i. e., we will use a superposition of two SU(2) states to solve the TDSE. Inspired by previous experience with Gaussian-based approaches to solve the TDSE for molecular Hamiltonians [30,31], as well as for spin-boson-type problems [32][33][34], and due to the entanglement entropy studies in [35] using two SU(2) states, we are confident that only a handful of suitable time-dependent basis states could be enough to achieve satisfactory agreement with exact quantum solutions if a full fledged variational approach is taken. In order to correctly account for more demanding quantum effects like self trapping, it will turn out that the multiplicity has to be increased, but it can still be kept below the total number of time-independent Fock states that has to be used in a full quantum calculation.
Furthermore, it is expected that the number of quantum trajectories needed for convergence will be much reduced as compared to semiclassical trajectory calculations that are based on multiple mean-field trajectories.
The presentation is structured as follows: In order to set the stage, in Sec. 2 we briefly review the mean-field approach, based on a single atomic coherent state (ACS), to the dynamics of the bosonic Josephson junction. At the end of this section, a special focus will be put on the stability analysis of the nonlinear classical phase space dynamics. In Sec.
3, we then choose an ansatz wave-function with non-trivial multiplicity, employing a small number of time-evolving atomic coherent states to represent the quantum beat dynamics (collapse and revival of the population imbalance) of the BH dimer. In a brief review of the quantum phase operator concept, we establish the relation between the phase difference in mean field and its quantum analog. This allows us to compare numerical results for phase space trajectories with the corresponding solution of the TDSE. We will cover a broad range of system parameters as well as initial conditions. It will turn out that there are cases, close to the equilibrium point of the classical dynamics, in which just two ACS will suffice to achieve reasonable agreement with exact results. Away from the classical equilibrium, the number of ACS will have to be increased, however. In the last section we give conclusions and an outlook on possible future work. Methodological details can be found in Section 5.

A. The Hamiltonian
The simplest Hamiltonian for the bosonic Josephson junction (two-site BH model) in normal ordered form readsĤ where the bosonic ladder operatorsâ j andâ † j with commutation relation [â j ,â † j ] =1 destroy, respectively create a particle (a bosonic atom) in the site labelled by the index j.
counts the number of particles per site andŜ =n 1 +n 2 is the total number operator and its expectation value S is a conserved quantity, becauseŜ commutes withĤ.
The (dimensionless) parameters U and J > 0 denote the strength of the on-site interaction, determined by the s-wave scattering length of the atomic species considered, and the tunneling amplitude, respectively. Later on, we will consider positive and negative values of U , corresponding to repulsive and attractive interaction between the atoms, respectively.

B. Mean-field dynamics
The evolution of the BH model is governed by the time-dependent Schrödinger equation for the wave-function |Ψ(t)⟩. Here as well as in the remainder of this paper we have set ℏ = 1. In order to solve for the dynamics, in the present section we will be following closely the mean-field work presented in [25].
We start the discussion, by recalling the eigenvalue equation of the annihilation operator in the formâ with the time-dependent Glauber coherent state [36] (displacement operator applied to the ground state) and time-dependent average particle number n j (t) and phase ϕ j (t) of the site indexed by j.
The position space representation of this state is a displaced Gaussian wavefunction [30].
The approximate mean-field dynamics can then be obtained by using an ansatz in terms of SU(2) coherent states (also refered to as atomic coherent states (ACS) [29]), defined by Here |0, 0⟩ is a shorthand notation for the direct product of two single-particle vacuum states and the time-dependent parameters and are the (normalized) population imbalance and the relative phase of the two sites, respectively [25]. The time-dependent particle number expectations at site j can take on fractional values. As a simple example, the reader may want to consider the case of S = 2 and initial z(0) = 1/2, for which n 1 (0) = 3/2 and n 2 (0) = 1/2.
If the system dynamics is governed by a harmonic oscillator Hamiltonian or a Rabi model (single harmonic mode coupled to a spin system), the use of the Glauber coherent states mentioned above is common [37,38]. In the present case, we opt for using the generalized coherent states (GCS) [39,40], which for two modes are the SU(2) coherent states introduced above, instead of a direct product of Glauber coherent states, however. This is because the former are better suited to describe particle number conserving dynamics, as the latter consist of a superposition of number states in the general case [41,42]. For Bose-Einstein condensates (BEC) this observation has also been made by Schachenmayer et al. [43], who showed that the multi-well Glauber coherent state ansatz is equivalent to the GCS ansatz only in the case of large particle numbers. Furthermore, it is worthwhile to note that the highly entangled GCS is the ground state of the "free-boson" model, i. e., the BH model with vanishing on-site interaction, U = 0 [7,[44][45][46].
where the entries of the vector ⃗ ξ are the complex parameters {ξ i }, which obey the "normal- The representation of the unit operator in these states has been used in [48] to establish an exact variational dynamics of the multi-mode Bose-Hubbard model. The number of independent real parameters of the GCS in the two-site case, M = 2, is three (two complex numbers minus the normalization condition mentioned above) but there is an overall phase factor that is irrelevant, however, so that we just remain with the two real parameters z and ϕ introduced above. In the case of arbitrary site numbers, the equations for the parameters ξ i are refered to as discrete nonlinear Schrödinger equation, which can be viewed as the discrete analog of the Gross-Pitaevskii equation for a BEC [6].
In the following, we will focus on the Josephson junction case. The mean-field equations for the real parameters z(t) and ϕ(t) are given by [25] which are equations of motion of non-rigid pendulum type [49][50][51]. A stationary solution of these coupled nonlinear equations is given by the equilibrium points (0, 2πn) with n = 0, ±1, ±2, . . . .
In the next step, we linearize the system of equations around one of the equilibrium points. The Jacobian matrix [52] at (z * , ϕ * ) = (0, 0) is given by and its eigenvalues are The so-called strength parameter is an appropriate parameter combination to be used frequently in the following. More details on the linearized mean-field equations around the stationary points can be found in A qualitative change in the mean field dynamics will occur, when the radicant in Eq. (13) changes sign, which happens at the critical value Λ SSB = −1, where the index SSB stands for spontaneous symmetry breaking [45]. If Λ > −1, both eigenvalues are imaginary, which indicates that the above equilibrium point is a stable one and the solution is symmetric around the origin, whereas Λ < −1 will lead to the emergence of another class of stable equilibrium points. The symmetry breaking solutions are located around the new stationary of the system of equations (10,11), where n ∈ Z [25]. The corresponding Jacobi matrix is given by and its eigenvalues are In [25] it is shown that the mean-field prediction for the onset of SSB based on single Glauber coherent states fails dramatically at small particle numbers. The mean-field result based on a single SU(2) coherent state does better than the Glauber state prediction at small S but is not exact. Both mean field predictions reproduce the full quantum result more faithfully at large values of S, however, as can be seen in Fig. 4 of [25].
A further conclusion that can be drawn from the mean-field equations (10,11) is the fact that the quantity is a constant of motion [25]. This leads to the existence of a parameter regime, in which the imbalance cannot become zero during an oscillation cycle and therefore, the average value of z will be nonzero. The condition for this macroscopic quantum self trapping (MQST) effect is E(z(0), ϕ(0)) > E(0, π) = JS. In terms of the strength parameter introduced above, the onset of self trapping is at [54] depending strongly on the initial position in phase space. In contrast to the case of SSB, the mean field MQST effect sets in at too large positive values of U (repulsive interaction), whereas mean field predicts the onset of SSB at too small values of |U | [25].

III. BEYOND MEAN FIELD DYNAMICS
Due to the shortcomings of the mean-field approach for U ̸ = 0, as there are the absence of collapses and revivals of the population imbalance [17], as well as failures in the prediction of the onset of MQST as well as SSB [25], we will now go beyond mean field by employing a multi-configuration ansatz for the solution of the TDSE.
We first give an explicit derivation of the equations of motion followed by a brief review of the phase operator concept, which is needed to display our quantum results. The parameter regimes of the results to be presented are put together in Table 1 We stress that all the parameters, compactly written as vectors A (with N entries) and ξ The corresponding Euler-Lagrange equations are given by where u k denotes one element of the set {A k , ξ k1 , ξ k2 } of 3N complex valued parameters in |Ψ⟩. For the coefficients, this leads to the equations of motion where For the coherent state parameters ξ jm (m = 1, 2), we get where and an analogous equation for the second index being 2.
Some numerical tricks to solve the highly non-linear, implicit equations of motion in Eqs. (23,25) have been devised in [55] for the case of Glauber coherent state basis functions.
Because of the restriction to M = 2 of the site number in the present investigation, at least for moderate particle numbers, the TDSE can also be solved easily by an expansion of the wave-function in (time-independent) Fock states, whose coefficients fulfill a (numerically more well-behaved) system of coupled linear first order differential equations. The number of Fock states that are required is determined by the particle number via S + 1. More details on the full quantum (Fock space) calculations, whose results will be refered to as exact quantum results in the following, are given in Appendix V B.
We stress that in all exact as well as beyond mean-field calculations to be presented, we take a single ACS as the initial condition of the dynamics. For the beyond mean-field calculations this means that a single element out of the set {A k } is nonzero initially, whereas all other elements will take nonzero values only in the course of time.

B. Brief Review of Phase Operator Concept
From the wave-function given by Eq. (20), we can calculate the time-dependent site populations by taking expectation values of the operators from (2) and from this the imbalance z between the two sites. For the analog of the relative phase ϕ, we use the quantum phase operator concept [9,56], leading to the expectation values of the cosine and sine of the phase operator and its sine square. In addition, the variance of the sine is defined by The normalization condition ⟨sin 2φ + cos 2φ ⟩ = 1 and the expectation of cos 2φ from [9] have been used to derive Eq. (29).
The relation between the sine of the classical phase variable, displayed in Fig. 1 and the expectation of the sine of the quantum phase is as can be derived by applying the operator in (28) to an ACS. For S → ∞, the prefactor on the RHS of the above equation becomes unity and the quantum and classical expressions become identical. Furthermore, in [9] it has been shown that the melting of coherence between the two sites is mirrored by the vanishing of the expectation of cosφ and the occurrence of large fluctuations of the corresponding variance.
In the following, we will focus on parameters on the border as well as inside of the most interesting regime, the so-called Josephson regime, for which the parameters fulfill the condition 1 < U S/(2J) ≪ S 2 [27]. Depending on the initial conditions, beyond mean-field effects can be observed in this case. In addition, we will also allow for negative values of the strength parameter Λ smaller than -1, in order to study the SSB case and will employ large positive Λ values close to the (mean-field) MQST regime.

C. Plasma oscillations
In the following, we first consider the case of small on-site interaction. In addition, also the initial imbalance shall first be small. In a second step, this imbalance shall be large at t = 0.

Small initial population imbalance
For small values of U , as well as of z, we only include two ACS in the Ansatz in Eq. (20), i. e., we use N = 2. Initially, ξ can be still be parameterized in analogy to the procedure of the previous section by (ξ 11 , ξ 12 ) = ( 1+z 1 2 , 1−z 1 2 e −iϕ 1 ) and (ξ 21 , ξ 22 ) = ( 1+z 2 2 , 1−z 2 2 e −iϕ 2 ). To highlight the changes that the inclusion of an additional basis state leads to, for the first SU(2) state, we use three different initial conditions, namely z 1,i ∈ {0.01, 0.05, 0.1}, ϕ 1 = 0, A 1 = 1. For the second SU(2) state, the initial values are identical and are fixed as The phase space trajectories in three different levels of approximation and for the three different initial conditions specified above are shown in Fig. 2, for S = 20 and an on-site interaction strength of U = 0.1J, implying U S/(2J) = 1. We stress that the expectation value of the sine of the phase operator is plotted on the y-axis. Its relation to the sine of the phase difference in the mean-field case is given in Eq. (31). In mean-field approximation all three initial conditions give rise to an ellipsoidal phase space pattern as shown in the previous section. Going beyond mean field by allowing for just one additional SU(2), we see a qualitatively different behaviour which corresponds to a beating of the population imbalance, here displayed by a spiraling motion that first moves inward and then outward for all three initial conditions. This is shown not to be an artefact by comparison to the full quantum solution, which shows an almost quantitative agreement with the ACS solution of multiplicity two. We stress that the choice of the initial phase of the second ACS is decisive for the quality of our beyond mean-field results. Choosing ϕ 2 to be zero, e.g., would lead to a spiraling in the wrong direction.
In the case of the smallest imbalance displayed in Fig. 2, the beating amplitude (the width of the blue ring) is smallest and the description of the quantum dynamics with a single classical (mean-field) trajectory is almost adequate, as it would be in the Rabi-oscillation regime, in which Λ ≪ 1, a case, we are not considering herein. As has been shown in [14], in order to cope with the collapse and revival of the population imbalance oscillations, a multitude of classical trajectories is needed, however. A ballpark number for the sample size in the Monte-Carlo integrations performed by Simon and Strunz is 10 4 . The truncated Wigner approximation based on a similar sampling procedure does not capture the revival oscillations, but a full fledged semiclassical approach is required to this end. To put our work in context, we stress, that to capture the quantum behaviour displayed in Fig. 2 almost quantitatively, we need only two "trajectories", i.e., two ACS. This dramatic reduction in basis size is due to the fact that in our present case, also the trajectories (the dynamical evolution of the basis function parameters) undergo the full variational procedure, i.e., they are not mean-field trajectories. We are thus loosing the intuitive appeal of a semiclassical method at the benefit of much smaller computational effort, although the calculation of the quantum trajectories is more involved than that of the mean-field ones. In passing, we note that in a comparison of the hierarchy of variational methods based on Glauber coherent states, applied to the anharmonic Morse potential in [30], the reduction in basis function size was counteracted by the (numerical) complexity of the solution of the variational equations of motion.

Large initial population imbalance
The small initial imbalance in the previous case has led to an incomplete collapse, i.e., the oscillation amplitude was still rather large in all three cases at all times, with only small beating amplitude. In order to suppress the total oscillation amplitude, i.e., to see very small oscillations at least temporarily, we have to allow for larger initial imbalances, which will be done next. For the case of U = 0.1J and with initial z = 0.5, the exact quantum dynamics for two different total particle numbers, S = 20 and S = 50, is shown in Fig.   3. Corresponding mean-field calculations (not shown) would display a closed single loop oscillation without any spiraling in (decrease of the oscillation amplitude). In the quantum case, however, we see an almost complete collapse of the amplitude, the larger the particle number. For the larger S, in addition, the population imbalance as well as the expectation of the sine of the phase operator stay around zero for a longer time (see also Fig. 4).
In Fig. 4 we display the time evolution of the sine of the phase operator and its variance as defined above. The results of the beyond mean-field approach and exact quantum calculations are compared. First, we observe the collapse and revival in the case of S = 20.
For S = 50, the maximum time considered is too short to observe the revival, however.  The occurrence of large amplitude oscillations in the variance has to be accounted for by an increase in the multiplicity, because in the single ACS case the relative phase is well-defined.
We note that both multiplicities are smaller than the total number of Fock states required, which is S + 1. Furthermore, the choice of the initial conditions for the initially unpopulated ACS is done in the random fashion explained in detail in [48].

D. Macroscopic quantum self trapping
In order to observe self trapping in the Josephson regime, i.e., the restriction of the population dynamics, such that the population on one side is always larger than on the other side the initial condition and/or the on-site interaction strength has to be changed.
From a mean-field argument the condition given in Eq. 19 has been derived, which is valid at all times. In the following, we will use z(0) = 0.5 and ϕ(0) = 0. This leads to Λ MQST ≈ 15.
We are choosing the total number of particles and the on-site interaction strength such that the actual value of Λ is just below the critical mean-field one, and that the classical dynamics therefore will not be trapped but the strength parameter is large enough for the quantum trajectory to be trapped at positive values of z [25]. For S = 20, we take U/J = 1.  In Fig. 5, the results for the phase space trajectories followed up to a total time of T = 50J are displayed. As dictated by our choice of parameters, the mean-field results do not display the MQST effect just yet. The quantum MQST has set in already, however. The fact that in the exact quantum results MQST happens for smaller coupling strengths than in mean field has also been reported in [25]. There it was found that the use of a single ACS does not allow one to observe this quantum effect (the reduction of the critical Λ value), however.
By consulting the red curves in Fig. 5, it can be seen that in order for the ACS-ansatz to show the correct quantum behavior (early onset of MQST) a non-trivial multiplicity has to be employed. In our present case this is N = 12 for the case S = 20 and N = 25 for the case S = 50. The high multiplicities needed are due to the fact that both the initial condition in z and also U are rather large. A single ACS will only give the exact quantum solution for U = 0, however.

E. Spontaneous symmetry breaking
So far, we have focused on positive on-site interaction strength. It was already observed on the mean-field level, however, that a spontaneous symmetry breaking is triggered by negative values of U beyond a certain threshold. The comparison of the mean-field with the full quantum solution as well as our multi-configuration ACS approach to this case will be in the focus of the present section.
Because the mean-field prediction for SSB is good for large particle numbers [25], in Fig.   6 we first consider the case of S = 50 and we take U/J = −0.12, leading to Λ < −1. Results of three different levels of approximation are again displayed: mean-field, ACS with small multiplicity (here N = 10) and full quantum. In the mean-field case, displayed in panel (a), we observe that the elliptic orbit for small deviations from the symmetry breaking equilibrium point (z SSB ≈ 0.94, ϕ SSB = 0), for larger displacements turns into a plectrum shaped orbit around the new stable fixed point (see also panel (b) of Fig. 1). As in the previous section, multi-configuration ACS with a small multiplicity of N = 10 displays the spiraling away from the mean-field orbit (the "quantum effect") in very faithful way. The further away from z SSB the initial condition is, the broader the range of the spiraling motion turns out to be, both in the ACS (panel (b)) and the exact results (panel (c)). The case of smaller particle number S = 20 and U/J = −0.15 leads to z SSB ≈ 0.71. The on-site interaction parameter lies just between the classically predicted onset of SSB and the quantum prediction. In the quantum case, it was shown that the SSB effect comes along with the switch from a unimodal to a bimodal distribution of the amplitudes in a Fock space expansion of the ground state of the Hamiltonian [45], which for large |U | becomes a so-called Schrödinger cat (NOON) state, in our notation a superposition proportional to |S, 0⟩+|0, S⟩. For smaller particle numbers, the onset of this effect, compared to the mean-field prediction, is pushed to larger absolute values of U (i. e., stronger attractive interaction), as shown in Fig. 4 of [25]. Thus for the parameters mentioned above, the phase space trajectories in the beyond mean-field case show a much different behavior than in the classical case. To visualize this behavior, we have calculated the Husimi transform [42] Q(z, ϕ) = ⟨Ω|Ψ⟩⟨Ψ|Ω⟩,  If one just wants to determine the transition from localized motion to delocalized motion beyond the mean-field prediction, it turned out above that only two ACS trajectories might be enough. We will show in the remainder of this section that in order to almost faithfully predict the occurence of the SSB transition in terms of |U | for different particle numbers, a multiplicity of N = 2 is indeed sufficient. To show this, in Fig to be sure that the motion is either confined to the right hand side of phase space (i.e., z > 0) or not and have used an interval nesting strategy to determine the onset of SSB.
These two results are then compared to the exact quantum ones (blue dots), calculated by monitoring the expansion of the ground state (GS) in terms of Fock states. If the magnitude of coefficients shows a bimodal structure, and ⟨GS| S 2 , S 2 ⟩ ≈ 0, the SSB range is reached [45]. The multi-configuration ACS calculations with N = 2 show a surpringly good agreement with the exact quantum results, even for small particle numbers. interactions, based on Gaussian wavepackets (Glauber coherent states) [58] as well as to the multi-configurational time-dependent Hartree-Fock method for bosons [59], although in the latter case, the employed basis functions are orthogonal. Furthermore, in contrast to the Glauber coherent states, the ACS used herein conserve particle number and are thus considered to be favorable in the present case [42]. In addition, we stress that in contrast to semiclassical methods that are based on Monte Carlo sampling of initial conditions for mean-field trajectories and require around 10 4 samples, here we can get satisfactory results with only a handful of variationally determined "trajectories". The semiclassical method employed by Tomsovic et al [8], requires an order of magnitude less mean-field trajectories (even in the 6 well case) than the initial value semiclassical method refered to above, but one has to find saddle points in a complexified phase space, which is a formidable task.
The parameter space that we have covered is characterized by the strength (and the sign) of the on-site interaction as well as by the total particle number and the initial population imbalance. Firstly, by taking into account one additional ACS, i. e., by employing a total of just two ACS, the beating of the population imbalance (as well as of the expectation of the sine of the phase operator) for small positive values of U can be reproduced almost quantitatively exactly, if the initial imbalance is rather small, i. e., if it is close to the classical equilibrium point at the origin of phase space. The choice of the initial phase variable of the second ACS was crucial to achieve this agreement. For larger initial imbalance, the number of ACS needed to achieve reasonable agreement with the exact quantum results has to be increased, with more and more states needed, the higher the total particle number.
Secondly, our focus was on the more demanding parameter regime of MQST. There, we could show that the use of more than ten ACS is necessary, if the quantum reduction compared to the mean-field value of the repulsive interaction strength at which MQST sets in is to be uncovered. As had been noticed before by Wimberger et al. [25], a single ACS is not enough to observe this effect. In the case of higher multiplicities N > 2, the choice of initial conditions for those ACS that are initially unpopulated (i., e., the ones, whose coefficients A k in Eq. (20) are zero) was done by the random sampling described in [48].
Lastly, for negative values of the on-site interaction and for large particle numbers, we observed a beating oscillation around the symmetry breaking equilibrium point, that still resembles the mean-field trajectory, with the only quantum effect being the spiraling in and out of the phase space trajectory. For small particle numbers, compared to the meanfield prediction, the symmetry breaking only occurs for larger attractive interaction in the quantum case, however [25]. The fact that the symmetry breaking is lost for parameters that would allow for symmetry breaking in mean-field theory is uncovered by using just two ACS. The new prediction of the onset of symmetry breaking in Fig. 8 is very close to the exact quantum result.
In future works, the fact that the addition of only a few generalized coherent state basis functions allows for the unraveling of quantum effects can be put to good use. A possible extension of the present work would be keeping the site number at two but allowing for more than just a single atomic species [60]. Furthermore, driven bosonic Josephson junctions show dynamical tunneling [61] and the addition of a decay term in one of the sites allows for a characteristic modulation of the self-trapping [62]. The description of these effects beyond for the population imbalance and the phase difference, valid around the phase space origin.
Employing the initial conditions z(0) = z 0 , ϕ(0) = 0, their solution is given by with the strength parameter Λ, defined in Eq. (14) of the main text.
If Λ > −1, we have the oscillatory solutions z(t) = z 0 cos(Ωt) (37) ϕ(t) = −z 0 √ 1 + Λ sin(Ωt) (38) with the plasma frequency We stress that for the specific choice of initial condition, the oscillation amplitude of ϕ depends on the strength parameter, while that of z does not.
If Λ < −1, we obtain which describes a solution like the ones displayed in panel (b) of Fig. 1, but only where the conditions |z|, |ϕ| ≪ 1 are still fulfilled. Away from that regime the hyperbolic solution is unphysical.

B. Exact quantum calculation
For the exact quantum results, we employ an expansion of the wave-function in terms of Fock states where the sum is taken over all the states {|F i ⟩} that emerge if a total of S particles is distributed over two sites. Due to the fact that one can place from zero up to S particles in, e.g., the first site, it is obvious that there are S + 1 different possibilities.
In order to completely specify the problem, the initial state has to be known, from which the b coefficients at t = 0 can be extracted. In the present work, we consider an initial state that is given in terms of a single ACS with parameters ξ 1 and ξ 2 . From the definition given in Eq. (9) taken for M = 2, by applying the binomial theorem, due to (â † i ) n |0⟩ = √ n!|n⟩, we find which is the Fock state expansion of the initial state, providing us with the sought for coefficients at t = 0.
The first option to evolve the wave-function over time would be to solve the coupled system of linear differential equations for the b coefficients that follows from the TDSE, e. g., by using a Runge-Kutta method or by matrix exponentiation (which in the present case of time-independent Hamiltonian turns out to be advantageous, because the matrix exponential has to be calculated only once, before the propagation loop is started). An alternative, second option, which is also numerically exact, would require diagonalising the BH Hamiltonian [16], e.g., in the Fock basis, see also Section 2.3.1 in [63]. The time evolution is then finally given by (ℏ = 1) where {E i } are the eigenenergies and the {|Φ i ⟩} are the eigenstates. The time-independent c-coefficients follow from the expansion of the initial wave-function in the eigenstates.
In both cases, the matrix elements of the Hamiltonian have to be set up. This does not pose a major challenge in case of small site numbers but in the general case it requires some clever way of creating and labeling of the Fock states, as described in a pedagogical way in [64].

CONFLICT OF INTEREST STATEMENT
were displayed here.