The simplest Hamiltonian for the bosonic Josephson junction (two-site BH model) in normal ordered form reads H ̂ = − J a ̂ 1 † a ̂ 2 + a ̂ 2 † a ̂ 1 + U 2 ∑ j = 1 2 a ̂ j † 2 a ̂ j 2 , (1) where the bosonic ladder operators a ̂ j and a ̂ j † with the commutation relation [ a ̂ j , a ̂ j † ] = 1 ̂ annihilate or create, respectively, a particle (a bosonic atom) in the site labeled by the index j. Furthermore, n ̂ j = a ̂ j † a ̂ j , j = 1,2 , (2) counts the number of particles per site and S ̂ = n ̂ 1 + n ̂ 2 is the total number-operator, and its expectation value S is a conserved quantity because S ̂ commutes with H ̂ .

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.

The evolution of the BH model is governed by the TDSE i | Ψ ̇ t 〉 = H ̂ | Ψ t 〉 , (3) 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 closely following the mean-field work presented in [25].

We start the discussion by recalling the eigenvalue equation of the annihilation operator in the form a ̂ j | α j t 〉 = α j t | α j t 〉 = n j t e i ϕ j t | α j t 〉 , (4) with the time-dependent Glauber coherent state [36] (displacement operator applied to the ground state) | α j t 〉 = e − | α j t | 2 / 2 e α j t a ̂ j † | 0 〉 (5) 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 wave-function [30].

The approximate mean-field dynamics can then be obtained by using an ansatz in terms of a single ACS [29], defined by | Ψ t 〉 = 1 S ! 1 + z t 2 a ̂ 1 † + 1 − z t 2 e − i ϕ t a ̂ 2 † S | 0,0 〉 = | S , 1 + z t 2 , 1 − z t 2 e − i ϕ t 〉 . (6) Here, |0, 0⟩ is a shorthand notation for the direct product of two single-particle vacuum states, and the time-dependent parameters z t = n 1 t − n 2 t S , (7) and ϕ t = ϕ 1 t − ϕ 2 t , (8) 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, we consider the case of S = 2 and initial z(0) = 1/2, for which n ₁(0) = 3/2 and n ₂(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 previously is common [37, 38]. In the present case, we opt for using the generalized coherent states (GCSs) [39, 40], which for two modes are the ACS introduced above, instead of a direct product of Glauber coherent states. 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 (BECs), 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–46].

The representation of the ACS in the last line in Eq. 6 is motivated by the general expression of a multimode generalized coherent state (total number of modes given by M) in the form [47] | S , ξ ⃗ 〉 = 1 S ! ∑ i = 1 M ξ i a i † S | 0,0 , … , 0 〉 , (9) where the entries of the vector ξ ⃗ are the complex parameters {ξ _i }, which obey the “normalization” condition ∑ i = 1 M | ξ i | 2 = 1 . The representation of the unit operator in these states has been used in [48]to establish an exact variational dynamics of the multi-mode BH model. The number of independent real parameters of the GCS in the two-site case, M = 2, is 3 (two complex numbers minus the normalization condition mentioned previously) but there is an overall phase factor that is irrelevant, however, so that we just remain with the two real parameters z and ϕ introduced previously. In the case of arbitrary site numbers, the equations for the parameters ξ _i are referred to as discrete non-linear Schrödinger equations, which can be viewed as the discrete analog of the Gross–Pitaevskii equation for a BEC [6].

In the remainder of this study, we will focus on the Josephson junction case. The mean-field equations for the real parameters z(t) and ϕ(t) are given by [25]: z ̇ = 2 J 1 − z 2 sin ⁡ ϕ ≔ f 1 , (10) ϕ ̇ = − 2 J z 1 − z 2 cos ⁡ ϕ − U S − 1 z ≔ f 2 , (11) which are equations of motion of non-rigid pendulum type [49–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 J = ∂ f 1 ∂ z z * , ϕ * ∂ f 1 ∂ ϕ z * , ϕ * ∂ f 2 ∂ z z * , ϕ * ∂ f 2 ∂ ϕ z * , ϕ * = 0 2 J − 2 J − S − 1 U 0 , (12) and its eigenvalues are λ ± = ± 2 J − 2 + U J − U S J . (13) The so-called strength parameter Λ = U S − 1 / 2 J , (14) is an appropriate parameter combination to be used frequently in the following section. More details on the linearized mean-field equations around the stationary points can be found in Section 5.1.

A qualitative change in the mean-field dynamics will occur when the radicand 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 aforementioned 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 point(s) z SSB , ϕ SSB = ± 1 − 1 Λ 2 , 2 π n , (15) of the system of Eqs 10, 11, where n ∈ Z [25]. The corresponding Jacobian matrix is J = 0 2 J Λ 2 − 2 J Λ 1 + Λ Λ 2 0 , (16) and its eigenvalues are λ ± = ± 2 J − Λ 4 − Λ Λ 2 Λ 2 . (17) If Λ < − 1, these are imaginary and the solution of the linearized equations around the SSB points is oscillatory. Both cases are displayed in Figure 1, with the left panel showing motion around the stable fixed point for U/J = 0.1 and the right panel showing the trajectories in case of U/J = −0.12, where the stable fixed point at the origin has turned into an unstable one and new stable fixed points appear at positive and negative values of z. For an experimental realization of this scenario, see [53]. In Section 5.1, an analytic expression for the linearized solution in panel B of Figure 1 at small values of z and ϕ is given.

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. However, both mean-field predictions reproduce the full quantum result more faithfully at large values of S, as can be seen in Figure 4 of [25].

It can be further concluded from the mean-field Eqs 10, 11 that the quantity E = U S 4 S − 1 z 2 − J S 1 − z 2 cos ⁡ ϕ , (18) is a constant of motion [25]. This leads to the existence of a parameter regime, in which the imbalance cannot become 0 during an oscillation cycle and, therefore, the average value of z will be non-zero. 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 previously, the onset of self-trapping is at [54] Λ MQST = 1 + 1 − z 2 0 cos ⁡ ϕ 0 z 2 0 / 2 , (19) 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 the mean-field theory predicts the onset of SSB at too small values of |U| [25].

As a step toward the exactness of the solution, we replace the wave-function of Eq. 6 by a linear combination of N time-dependent SU(2) coherent states, written as in the general SU(M) case of Eq. 9, leading to | Ψ t 〉 = ∑ k = 1 N A k t | S , ξ k 1 t , ξ k 2 t 〉 . (20) We stress that all the parameters, compactly written as vectors A (with N entries) and ξ (with 2N entries), are time-dependent and complex-valued. Their (non-linear) equations of motion, again derived from the TDVP, in the general case of arbitrary multiplicity N as well as the site number M, have been given in matrix form in the appendix of [48].

For being self-contained, here, we explicitly review the variational procedure for the double-well problem. With the trial state from Eq. 20 the Lagrangian L ≔ i ⟨ Ψ | ∂ t | Ψ ⟩ − ⟨ Ψ | H ̂ | Ψ ⟩ takes the explicit form L = ∑ k , j = 1 N A k * A ̇ j 〈 ξ ⃗ k | ξ ⃗ j 〉 + i S ∑ k , j = 1 N A k * A j ∑ i = 1 2 ξ k * ξ ̇ j i 〈 ξ k ′ ⃗ | ξ j ′ ⃗ 〉 − ∑ k , j = 1 N A k * A j − J S ξ k 1 * ξ j 2 + ξ k 2 * ξ j 1 〈 ξ k ′ ⃗ | ξ j ′ ⃗ 〉 + U 2 S S − 1 ξ k 1 * 2 ξ j 1 2 + ξ k 2 * 2 ξ j 2 2 〈 ξ k ′ ′ ⃗ | ξ j ′ ′ ⃗ 〉 . (21) The corresponding Euler-Lagrange equations are given by ∂ L ∂ u k * − d d t ∂ L ∂ u ̇ k * = 0 , (22) 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 i ∑ j = 1 N A ̇ j 〈 ξ k ⃗ | ξ j ⃗ 〉 + i S ∑ j = 1 N A j ∑ i = 1 2 ξ k i * ξ ̇ j i 〈 ξ k ′ ⃗ | ξ j ′ ⃗ 〉 − ∂ H ∂ A k * = 0 , (23) where ∂ H ∂ A k * = ∑ j = 1 N A j − J S ξ k 1 * ξ j 2 + ξ k 2 * ξ j 1 〈 ξ k ′ ⃗ | ξ j ′ ⃗ 〉 + U 2 S S − 1 ∑ i = 1 2 ξ k i * 2 ξ j i 2 〈 ξ k ′ ′ ⃗ | ξ j ′ ′ ⃗ 〉 . (24) For the coherent state parameters ξ _jm (m = 1, 2), we get i S ∑ j = 1 N A k * A ̇ j ξ j m ⟨ ξ k ′ ⃗ | ξ j ′ ⃗ ⟩ + ∑ j = 1 N A k * A j ξ ̇ j m ⟨ ξ k ′ ⃗ | ξ j ′ ⃗ ⟩ + S − 1 ∑ j = 1 N A k * A j ∑ i = 1 2 ξ k i * ξ ̇ j i ξ j m ⟨ ξ k ′ ′ ⃗ | ξ j ′ ′ ⃗ ⟩ − ∂ H ∂ ξ k m * = 0 , (25) where ∂ H ∂ ξ k 1 * = ∑ j = 1 N A k * A j − J S ξ j 2 〈 ξ k ′ ⃗ | ξ j ′ ⃗ 〉 − J S S − 1 ξ k 1 * ξ j 2 + ξ k 2 * ξ j 1 ξ j 1 〈 ξ k ′ ′ ⃗ | ξ j ′ ′ ⃗ 〉 + U S S − 1 ξ k 1 * ξ j 1 2 〈 ξ k ′ ′ ⃗ | ξ j ′ ′ ⃗ 〉 + U 2 S S − 1 S − 2 ∑ i = 1 2 ξ k i * 2 ξ j i 2 ξ j 1 〈 ξ k ′ ′ ′ ⃗ | ξ j ′ ′ ′ ⃗ 〉 , (26) 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 basis state 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 required is determined by the particle number via S + 1. More details on the full quantum (Fock space) calculations, whose results will be referred to as exact quantum results, are given in Section 5.2.

We stress that in all exact and 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 non-zero initially, whereas all other elements will take non-zero values only in the course of time.

From the wave-function given by Eq. 20, we can calculate the time-dependent site populations by taking expectation values of the operators from Eq. 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 ⟨ cos ϕ ̂ ⟩ = ⟨ a ̂ 2 † a ̂ 1 + a ̂ 2 a ̂ 1 † ⟩ 2 2 n ̂ 1 n ̂ 2 + n ̂ 1 + n ̂ 2 , (27) ⟨ sin ϕ ̂ ⟩ = i ⟨ a 1 † a 2 − a 2 † a 1 ⟩ 2 2 n ̂ 1 n ̂ 2 + n ̂ 1 + n ̂ 2 , (28) ⟨ sin 2 ϕ ̂ ⟩ = 1 2 − ⟨ a ̂ 2 † a ̂ 1 2 + a ̂ 2 a ̂ 1 † 2 ⟩ 2 2 n ̂ 1 n ̂ 2 + n ̂ 1 + n ̂ 2 , (29) of the cosine and sine of the phase operator and its sine square. In addition, the variance of the sine is defined by Δ sin ϕ ̂ ≔ 〈 sin 2 ϕ ̂ 〉 − 〈 sin ϕ ̂ 〉 2 . (30) 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 Figure 1, and the expectation of the sine of the quantum phase is ⟨ sin ϕ ̂ ⟩ = S 1 − z 2 S S − 1 1 − z 2 + 2 S sin ⁡ ϕ , (31) as can be derived by applying the operator in Eq. 28 to an ACS. For S → ∞, the prefactor on the RHS of the aforementioned 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 section, we will focus on parameters on the border and inside of the most interesting regime, the so-called Josephson regime [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 use large positive Λ values close to the (mean-field) MQST regime.

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

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 section, we will use z(0) = 0.5 and ϕ(0) = 0. This leads to Λ_MQST ≈ 15. We choose 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.2 and for S = 50, we take U/J = 0.53, leading to Λ ≈ 14.4 and Λ ≈ 13.0, respectively. Both values are deep inside the Josephson regime.

In Figure 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. However, the quantum MQST has set in already. The fact that in the exact quantum results, MQST happens for smaller coupling strengths than in mean-field has also been reported in [25]. 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). By observing the red curves in Figure 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 (displayed in panel A) and N = 25 for the case S = 50 (displayed in panel B). 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 result for U = 0.

So far, we have focused on the case of positive on-site interaction strength. However, it was already observed on the mean-field level 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 and our multi-configuration ACS approach for this case will be the focus of the present section.

Because the mean-field prediction for SSB is good for large particle numbers [25], in Figure 6, we first consider the case of S = 50 and we take U/J = −0.12, leading to Λ < − 1. The 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 Figure 1). As in the previous section, multiconfiguration ACS with a small multiplicity of N = 10 displays the spiraling away from the mean-field orbit (the “quantum effect”) in a 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 numbers 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 of the amplitude from a unimodal to a bimodal distribution 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 Figure 4 of [25]. Thus, for the parameters mentioned previously, 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 , ϕ = 〈 Ω | Ψ 〉 〈 Ψ | Ω 〉 , (32) with |Ω⟩ = |S, z, ϕ⟩. This function is localized if the time-evolved quantum state is localized around the stable fixed point, whereas it is delocalized otherwise. Taking snapshots of its dynamics for U/J = −0.15, displayed in panels A to C of Figure 7, a delocalization of the dynamics can be observed, which is very different from the mean-field prediction, which is already deep in the SSB regime. Using just two ACSs, we can thus unravel this quantum effect dynamically, without having to calculate the exact ground state. Increasing the absolute value of the onsite interaction to |U|/J = 0.19, according to Eq. 5, the new stable fixed points are moving towards larger |z| and we are choosing an initial condition with a small displacement away from the one with ϕ = 0 and a positive value of z. The beyond mean-field result (with N = 2) now also shows a restriction of the dynamics to positive values of z, i.e., it displays the phenomenon of SSB. This fact can be observed in panels D to F of Figure 7. The fact that there is no motion from right to left, when |U| is large, is due to the large energy barrier that has to be overcome in order to go from z ≈ 1, i.e., approximately the state |S, 0⟩ to z ≈ − 1, i.e., approximately the state |0, S⟩ [57].

If one just wants to determine the transition from delocalized to localized motion beyond the mean-field prediction, it turned out 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 occurrence of the SSB transition in terms of |U| for different particle numbers, a multiplicity of N = 2 is indeed sufficient. To show this, in Figure 8, we display the onset of SSB for different values of S, as predicted by mean-field (using a single ACS), yielding the hyperbolic dependence |U|/2J = 1/(S − 1) depicted by the solid line, to the result for this onset from a multi-configurational calculation with N = 2 (red diamonds). To determine the location of the red diamonds, we have propagated the dynamics up to large enough times (Jt = 1000) 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 those of 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 G S | S / 2 , S / 2 ≈ 0 , the SSB range is reached [45]. The multi-configuration ACS calculations with N = 2 show a surprisingly good agreement with the exact quantum results, even for small particle numbers.

From the Jacobi matrix in Eq. 12 we read off the linearized equations of motion z ̇ = 2 J ϕ , (33) ϕ ̇ = − 2 J + U S − 1 z (34) for the population imbalance and the phase difference, valid around the phase space origin. Employing the initial conditions z(0) = z ₀, ϕ(0) = 0, their solution is given by z t = z 0 2 e 2 J − 1 + Λ t + e − 2 J − 1 + Λ t , (35) ϕ t = z 0 2 − 1 + Λ e 2 J − 1 + Λ t − e − 2 J − 1 + Λ t , (36) 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 Ω = 2 J 1 + Λ . (39) 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 z t = z 0 ⁡ cosh 2 J − 1 + Λ t , (40) ϕ t = z 0 − 1 + Λ sinh 2 J − 1 + Λ t , (41) which describes a solution like the ones displayed in panel B of Figure 1, but only where the conditions |z|, |ϕ|≪ 1 are still fulfilled. Away from that regime the hyperbolic solution is unphysical.

For the exact quantum results, we employ an expansion of the wave-function in terms of Fock states | Ψ t 〉 = ∑ i = 0 S b i t | F i 〉 , (42) 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 ξ ₁ and ξ ₂. From the definition given in Eq. 9 taken for M = 2, by applying the binomial theorem, due to ( a ̂ i † ) n | 0 〉 = n ! | n 〉 , we find | Ψ 0 〉 = 1 S ! ξ 1 a ̂ 1 † + ξ 2 a ̂ 2 † S | 0,0 〉 = ∑ j = 0 S S ! S − j ! j ! ξ 1 S − j ξ 2 j | S − j , j 〉 , (43) 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 i b ̇ j t = ∑ i = 0 S ⟨ F j | H ̂ | F i ⟩ b i t , (44) 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) | Ψ t 〉 = ∑ i = 0 S c i ⁡ exp − i E i t | Φ i 〉 , (45) 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].