# Exact variational dynamics of the multi-mode Bose-Hubbard model based on
SU($M$) coherent states

Yulong Qiao  Max-Planck-Institut für Physik Komplexer Systeme, Nöthnitzer Str.
38, D-01187 Dresden, Germany  Institut für Theoretische Physik, Technische
Universität Dresden, D-01062 Dresden, Germany  Frank Grossmann  Institut für
Theoretische Physik, Technische Universität Dresden, D-01062 Dresden, Germany

###### Abstract

We propose a variational approach to the dynamics of the Bose-Hubbard model
beyond the mean field approximation. To develop a numerical scheme, we use a
discrete overcomplete set of Glauber coherent states and its connection to the
generalized coherent states studied in depth by Perelomov [A. Perelomov,
Generalized Coherent States and Their Applications, Springer-Verlag (Berlin,
1986)]. The variational equations of motion of the generalized coherent state
parameters as well as of the coefficients in an expansion of the wavefunction
in terms of those states are derived and solved for many-particle problems
with large particle numbers $S$ and increasing mode number $M$. For $M=6$ it
is revealed that the number of complex-valued parameters that have to be
propagated is more than one order of magnitude less than in an expansion in
terms of Fock states.

##  I Introduction

The mathematical foundation of the phase space formulation of physical systems
with Lie group symmetries has been considerably widened by the works of Brif
and Mann Brif and Mann (1998, 1999). Based on their progress, the description
of the Bose-Hubbard (BH) dynamics in terms of phase space distributions has
received new impetus from the works of Korsch and collaborators Trimborn  _et
al._ (2008, 2009). Therein, equations of motion for the $P$\- as well as the
$Q$-function of quantum optics Zhang  _et al._ (1990); Gardiner and Zoller
(2004) have been derived and solved for small site numbers $M$. The basis
functions that were found appropriate for the treatment of the particle number
conserving dynamics are so-called SU($M$) coherent states (CS), also referred
to as generalized coherent state (GCS) Perelomov (1986). These have been
investigated and favored in the same context by Buonsante and Penna in their
review Buonsante and Penna (2008), whose focus is on variational mean field
methods. A more recent review with a focus on SU(2) CS, introduced as atomic
coherent states in Arecchi  _et al._ (1972), is given in Klimov  _et al._
(2017). An extension of the formalism towards dissipative Lindblad type
equations in terms of $P$ functions has been given in Merkel  _et al._ (2020).

In parallel, the use of discrete, complete von Neumann type sets of the more
“standard” Glauber coherent states, whose position representation are Gaussian
wavepackets, has been applied to a flurry of different physical situations,
ranging from electron dynamics in atoms to nuclear dynamics in molecules as
well as to nonadiabatic (combined electron-nuclear motion), as reviewed in
Grossmann (1999); Miller (2001); Thoss and Wang (2004); Kay (2005); Richings
_et al._ (2015); Werther  _et al._ (2021). The equations of motion of the
coherent state parameters as well as of the coefficients in the expansion of
the wavefunction in terms of those states are usually derived from a
variational principle and possibly undergo additional approximations. It has
turned out in numerical investigations that the use of a surprisingly small
number of CS basis functions leads to converged results, e.g., in spin-boson
type problems tackled by the so-called multi Davydov D2 Ansatz Wang  _et al._
(2017); Werther and Grossmann (2020a, b), or in the method called coupled
coherent states for indistinguishable bosons Green and Shalashilin (2019).

Furthermore, in Lee and Fischer (2014), a generalization of the multi-
configuration time-dependent Hartree method for bosons (MCTDHB) Alon  _et al._
(2008) based on McLachlan’s variational principle has been given. The time-
dependent permanents used there are based on orthogonal orbitals, however.
Previous approaches that have used Glauber coherent state based semiclassical
propagators for BH dynamics have been restricted to small mode numbers Simon
and Strunz (2014); Ray  _et al._ (2016), as is also the case for semiclassical
approaches based on SU($M$) CS Viscondi and de Aguiar (2011); Khripkov  _et
al._ (2013), whereas mean field approaches as the ones discussed in Buonsante
and Penna (2008) rely just on a single basis function (i.e. trivial
multiplicity).

It has been worked out in Rau  _et al._ (2010a, b), that Gaussian CS are
promising basis functions also for full fledged dynamical calculations for BH
dynamics beyond the semiclassical propagation or the Gross-Pitaevskii level
Polkovnikov  _et al._ (2002). Similarly, imaginary time propagation to find
the ground state as well as real-time propagation to extract the excitation
spectrum using Gaussians has been performed recently for BH systems Guaita
_et al._ (2019). In addition, in the two mode case, the use of a single time-
dependent atomic coherent state has been promoted in Wimberger  _et al._
(2021). These successes as well as that of the CS basis functions alluded to
above, leads us to investigate the question if also for the fixed particle
number generalized coherent states multi-mode BH dynamics can be treated by a
numerically exact variational approach, based on a multitude of those GCS.

The presentation is organized as follows: First, in Section II, we review the
one-dimensional BH model and the relation between the Glauber (field) coherent
states Louisell (1990) and the (generalized) SU($M$) coherent states Perelomov
(1986). This latter review is necessary in order to get a handle on the
discretization of the representation of unity in terms of these states. In
Section III, we then derive the variational equations of motion for the GCS
parameters using the Lagrangian of the time-dependent variational principle.
As a proof of principle in Sect. IV we finally solve the dynamical problem for
several realizations of the BH model with different (relatively large) mode
number $M$ and particle number $S>M$ and compare to results gained from an
expansion in Fock states where this is still feasible. Conclusions are given
in Sect. V, while in the appendix some properties of the GCS that are needed
in the main text, as well as the matrix form of the variational equations and
a detailed convergence study are gathered.

##  II Bose-Hubbard model and completeness of SU($M$) coherent states

In the following we are constructing a discrete set of GCS for use in a
variational approach to the dynamics of the BH model, which is going to be
reviewed first.

###  II.1 Bose-Hubbard model in one dimension

In this presentation, we are aiming at a treatment in terms of the variational
principle of the dynamics of the one-dimensional linear chain BH model Jaksch
and Zoller (2005)

$\displaystyle\hat{H}=-J\sum_{j=1}^{M-1}\left[(\hat{a}_{j}^{\dagger}\hat{a}_{j+1}+h.c.)\right]+\frac{U}{2}\sum_{j=1}^{M}\hat{a}_{j}^{{\dagger}2}\hat{a}_{j}^{2}+\frac{K}{2}\sum_{j=1}^{M}(j-j_{0})^{2}\hat{a}_{j}^{\dagger}\hat{a}_{j},$
(1)

written here in normal-ordered form in terms of bosonic creation and
annihilation operators, $\hat{a}_{j}^{\dagger}$ and $\hat{a}_{j}$, which
fulfill the commutation relation
$[\hat{a}_{i},\hat{a}_{j}^{\dagger}]=\hat{\delta}_{ij}$ and their action on
number states $\\{|n_{j}\rangle\\},n_{j}=0,1,\dots$ is given by

$\hat{a}_{j}|n_{j}\rangle=\sqrt{n_{j}}|n_{j}-1\rangle,\\\
\hat{a}_{j}^{\dagger}|n_{j}\rangle=\sqrt{n_{j}+1}|n_{j}+1\rangle.$ (2)

The hopping matrix element $J$ is determined by the tunneling probability
between nearest neighbors, whereas $U$ denotes the on-site interaction, and
$K$ is due to an external, harmonic trapping potential (centered around
$j_{0}$), which will, however, be set to zero for all the numerical results to
be presented below. The above Hamiltonian is a widely studied model for the
dynamics of spinless particles in optical lattices Bloch  _et al._ (2008). An
experimental realization of a quantum phase transition has been reached by
varying the model’s parameters Greiner  _et al._ (2002). Furthermore, Kolovsky
has given a review of the spectral properties of the BH Hamiltonian in the
light of quantum chaology Kolovsky (2016).

In Trimborn  _et al._ (2008) it has been argued that an analysis in terms of
flat space and the use of the corresponding Glauber coherent states to
describe particle number preserving BH dynamics is inadequate. The reason
given, in the case of $M$ sites (or modes), is that the dynamical group is
spanned by the normally ordered operators

$\hat{a}_{i}^{\dagger}\hat{a}_{j},\qquad i,j\in\\{1,2,\dots M\\}$ (3)

and is equivalent to the special unitary group SU($M$). Therefore the
corresponding GCS should be employed to investigate the dynamics. We note in
passing that, in terms of the normally ordered operators from Eq. (3), the
operators $\hat{a}_{j}^{{\dagger}2}\hat{a}_{j}^{2}$ from Eq. (1) are written
as $\hat{n}_{j}(\hat{n}_{j}-1)$, with the number operators
$\hat{n}_{j}=\hat{a}_{j}^{\dagger}\hat{a}_{j}$.

###  II.2 Discrete grid of generalized coherent states

An insightful didactic discussion of the group theoretical construction of
coherent states is given in Zhang  _et al._ (1990). Starting by revisiting the
simple case of Glauber field coherent states as the coherent states of the
Heisenberg Weyl group, more complex situations are discussed, where the group
theoretical approach will be more advantageous for the understanding of the
topological structure of the corresponding coherent states than in the simple
case. An open question remains as to the completeness of a discrete set of
coherent states in the general case, however.

The Glauber coherent states $\left|{z}\right\rangle$ (single site case) can be
defined via the eigenvalue equation

$\displaystyle\hat{a}\left|{z}\right\rangle$
$\displaystyle=z\left|{z}\right\rangle$ (4)

with the complex eigenvalue

$\displaystyle z$ $\displaystyle=\frac{\gamma^{1/2}q+{\rm
i}\gamma^{-1/2}p}{\sqrt{2}}.$ (5)

in terms of position and momentum (phase space) and width parameter $\gamma$.
These states consist of a Poissonian superposition of number states Louisell
(1990).

Exactly fifty years ago, two independent contributions have proven the
statement that a discretized version of the unit operator is given by

$\sum_{k,l}|z_{k}\rangle({\mathbf{\Omega}}^{-1})_{kl}\langle z_{l}|=\hat{I},$
(6)

with the overlap matrix ${\mathbf{\Omega}}$ with elements

$\Omega_{kl}=\langle z_{k}|z_{l}\rangle,$ (7)

if the requirement

$z_{k}=\beta(m+{\rm i}n),$ (8)

with $k=(m,n)$ where $m,n$ are integers and $0<\beta\leq\sqrt{\pi}$, for the
spacing of the grid points is met Perelomov (1971); Bargmann  _et al._ (1971).
Physically, this means that the cells in $(p,q)$ phase space that are spanned
by the grid points must have an area less than or equal to the Planck cell
area of $2\pi$ in the present units (where $\hbar=1$). The limiting case of
$\sqrt{\pi}$ for the spacing of the grid points has already been postulated
long before by von Neumann von Neumann (1955); Boon and Zak (1978).

As mentioned above, for the BH model, generalized SU($M$) coherent states are
the appropriate coherent states, staying within the particle number conserving
subspace of the dynamics. In order to get a handle on their completeness, we
use the expansion of a Glauber coherent state in terms of the generalized CS.
The representation of the multimode GCS that is most appropriate to this end
is Buonsante and Penna (2008)

$\displaystyle|S,\vec{\xi}\rangle=\frac{1}{\sqrt{S!}}\Big{(}\sum_{i=1}^{M}\xi_{i}a_{i}^{\dagger}\Big{)}^{S}|0,0,\cdots,0\rangle,$
(9)

where $|0,0,\cdots,0\rangle$ denotes the multi-mode vacuum state and $S$ is
the number of bosons of the GCS. Some properties of the GCS that will be
needed later-on are gathered in Appendix A. The set of complex numbers
$\\{\xi_{i}\\}$ are characteristic parameters of the GCS, and they satisfy the
normalization condition

$\displaystyle\sum^{M}_{i=1}|\xi_{i}|^{2}=1,$ (10)

where $M$ represents the number of different modes present in the BH model of
Eq. (1). Together with the total number of bosons $S$ (not restricted to, but
herein chosen to be larger than $M)$, this determines the dimension of the
Hilbert space

$\displaystyle C^{F}={M+S-1\choose M-1}={M+S-1\choose
S}=\frac{(M+S-1)!}{S!(M-1)!}$ (11)

spanned by the Fock states Kolovsky (2016). The relation of the above
definition of the GCS to the group theoretical formulation for $M>2$ is
obscured by the necessity to disentangle exponentiated operators Buonsante and
Penna (2008). It is helpful to write out the definition in Eq. (9) explicitly
for small numbers $M=2<S$ to convince oneself of the fact that the GCS is a
superposition of Fock states with coefficients to be determined from a
generalized binomial formula.

As shown in Buonsante and Penna (2008), and as can be proven by Taylor
expansion of the exponential function, the multimode Glauber coherent state
$|Z\rangle=\prod_{i=1}^{M}|z_{i}\rangle$, and the GCS are related by

$\displaystyle|Z\rangle\ $
$\displaystyle=e^{-\frac{1}{2}\sum_{i=1}^{M}|z_{i}|^{2}}e^{\sum_{i=1}^{M}z_{i}a_{i}^{\dagger}}|0,0,\cdots,0\rangle$
$\displaystyle=e^{-\frac{1}{2}\sum_{i=1}^{M}|z_{i}|^{2}}\sum_{S=0}^{\infty}\frac{1}{S!}(\sum_{i=1}^{M}z_{i}a_{i}^{\dagger})^{S}|0,0,\cdots,0\rangle$
$\displaystyle=e^{-\frac{\tilde{N}}{2}}\sum_{S=0}^{\infty}\frac{\tilde{N}^{\frac{S}{2}}}{\sqrt{S!}}\frac{1}{\sqrt{S!}}\Big{(}\sum_{i=1}^{M}{\frac{z_{i}}{\sqrt{\tilde{N}}}}a_{i}^{\dagger}\Big{)}^{S}|0,0,\cdots,0\rangle$
$\displaystyle=e^{-\frac{\tilde{N}}{2}}\sum_{S=0}^{\infty}\frac{\tilde{N}^{\frac{S}{2}}}{\sqrt{S!}}|S,\vec{\xi}\rangle.$
(12)

Here, $\tilde{N}=\sum_{i=1}^{M}|z_{i}|^{2}$ denotes the average number of
bosons in $|Z\rangle$. We note that the relationship between $\xi_{i}$ and
$z_{i}$ is $\xi_{i}=\frac{z_{i}}{\sqrt{\tilde{N}}}$. Thus it is natural to
construct a one-to-one map between the sets $\\{\xi_{k,i}\\}$ and
$\\{z_{k,i}\\}$ (here the first index denotes the basis function
discretization index and second index denotes the mode index) to be used in a
discretized version of the unit operator.

To this end, we first generalize the completeness relation for the single site
case from Eq. (6) for multimode Glauber coherent states: If the complex grid
$L_{P,Q}^{M}$ which lies in the multidimensional complex $Z$-plane satisfies
the multidimensional Planck cell condition, the multimode coherent states will
form an overcomplete set and they obey the closure relation

$\displaystyle\sum_{k,l}A_{k,l}|Z_{k}\rangle\langle Z_{l}|=\hat{I},\indent
Z_{k},Z_{l}\in L_{P,Q}^{M},$ (13)

where

$\displaystyle\bm{A}=\bm{\Omega}^{-1},\Omega_{kl}=\langle Z_{k}|Z_{l}\rangle,$
(14)

with the multi-dimensional direct product CS
$|Z_{k}\rangle=|z_{k1},z_{k2},\cdots z_{kM}\rangle$. A heuristic proof of the
closure relation can be obtained in the following way: Given any two coherent
states $|Z_{m}\rangle$ and $|Z_{n}\rangle$ and inserting unity in the multi
degree of freedom case yields

$\displaystyle\langle Z_{m}|Z_{n}\rangle$
$\displaystyle=\sum_{k,l=1}^{N}A_{kl}\langle Z_{m}|Z_{k}\rangle\langle
Z_{l}|Z_{n}\rangle$ $\displaystyle=\sum_{k,l=1}^{N}\langle Z_{m}|Z_{k}\rangle
A_{kl}\langle Z_{l}|Z_{n}\rangle$ (15)

This equation has the matrix form $\bm{\Omega}=\bm{\Omega}\bm{A}\bm{\Omega}$
with the overlap matrix ${\bf\Omega}$ from Eq. (14) above. With the definition
of the matrix ${\bf A}$ as the inverse of the overlap matrix, the matrix
equation becomes an identity.

From this, we infer that if a set of states is complete in the whole Hilbert
space, these states are also complete for any subspaces which belong to the
Hilbert space. In our present case, we take the subspace to be the one
consisting of multimode Fock states
$\\{|n_{S}\rangle=|n_{1},n_{2},\cdots\rangle|n_{1},n_{2},\cdots\in\mathcal{N},\langle
n_{S}|\hat{N}|n_{S}\rangle=S,\hat{N}=\sum_{i}\hat{n}_{i}\\}$ with fixed number
of bosons $S$.

Since the Glauber coherent state can be written as a sum of SU($M$) coherent
states, according to

$\displaystyle|z_{k1},z_{k2},\cdots,z_{kj},\cdots\rangle=e^{-\frac{\tilde{N}}{2}}\sum_{S=0}^{\infty}\frac{\tilde{N}^{\frac{S}{2}}}{\sqrt{S!}}|S,\xi_{k1},\xi_{k2},\cdots,\xi_{kj},\cdots\rangle,$
(16)

there is a map
$z_{kj}\rightarrow\xi_{kj}=\frac{z_{kj}}{\sqrt{\sum_{i=1}^{M}|z_{ki}|^{2}}}$.
Unlike the $|z_{k1},z_{k2},\cdots,z_{kj},\cdots\rangle$ where the parameters
$z_{k1},z_{k2},\cdots,z_{kj},\cdots$ are mutually independent, the
$\xi_{k1},\xi_{k2},\cdots,\xi_{kj},\cdots$ are correlated due to the
normalization factor $\frac{1}{\sqrt{\sum_{i=1}^{M}|z_{ki}|^{2}}}$, however.
It is obvious that if the $z_{k1},z_{k2},\cdots,z_{kj},\cdots$ form a complex
grid fulfilling the completeness criterion, then the corresponding
$\xi_{k1},\xi_{k2},\cdots,\xi_{kj},\cdots$ will also form a complete basis

$\displaystyle\Big{\\{}|S,\xi_{k1},\xi_{k2},\cdots,\xi_{kj},\cdots\rangle|\xi_{kj}=\frac{z_{kj}}{\sqrt{\sum_{i=1}^{M}|z_{ki}|^{2}}},z_{kj}\in
L_{P,Q}^{M}\Big{\\}}.$ (17)

for the subspace
$\\{|n_{S}\rangle=|n_{1},n_{2},\cdots\rangle|n_{1},n_{2},\cdots\in\mathcal{N},\langle
n_{S}|\hat{N}|n_{S}\rangle=S,\hat{N}=\sum_{i}\hat{n}_{i}\\}$.

We now consider the two-mode case as an instructive example. The two initially
independent complex grids for $z_{k1}$ and $z_{k2}$ are shown in Fig. 1.
Taking samples from these two grids for $z_{k1}$ and $z_{k2}$ randomly at the
same time, every produced pair of $\\{z_{k1},z_{k2}\\}$ forms a two-mode
Glauber coherent state $|z_{k1},z_{k2}\rangle$. The $\\{z_{k1},z_{k2}\\}$
correspond to a set of couples of parameters
$\\{\xi_{k1},\xi_{k2}\\}=\\{\frac{z_{k1}}{\sqrt{|z_{k1}|^{2}+|z_{k2}|^{2}}},\frac{z_{k2}}{\sqrt{|z_{k1}|^{2}+|z_{k2}|^{2}}}\\}$
for the SU(2) coherent states, but due to the mutual correlation between
$\xi_{k1}$ and $\xi_{k2}$, they cannot be presented on two independent complex
planes. We are now rewriting the GCS parameters in the form

$\displaystyle\xi_{k1}$ $\displaystyle=\cos\left(\frac{\theta_{k}}{2}\right)$
(18) $\displaystyle\xi_{k2}$
$\displaystyle=\sin\left(\frac{\theta_{k}}{2}\right)e^{i\phi_{k}},$ (19)

where $0\leq\theta_{k}\leq\pi$ and $0\leq\phi_{k}\leq 2\pi$ are the angles and
the relative phases, respectively. Because we have just two angles for every
$k$ and because of the factor of 1/2 that comes along with the angle
$\theta_{k}$, this allows us to represent the GCS parameters as points on the
northern hemisphere of the Bloch sphere.

Figure 1: The complex grids for $z_{1}$ and $z_{2}$ in the case of two modes.
They are centered around $\\{z_{01},z_{02}\\}=\\{\xi_{01},\xi_{02}\\}$ given
by the GCS parameters of the initial state.

For a visualization, we first generated 50 random pairs of
$\\{z_{k1},z_{k2}\\}$. Using the angles $\\{\theta_{k},\phi_{k}\\}$, the
corresponding 50 pairs of $\\{\xi_{k1},\xi_{k2}\\}$ can be expressed via
$\\{\cos(\theta_{k}/2)\cos(\phi_{k}),\cos(\theta_{k}/2)\sin(\phi_{k}),\sin(\theta_{k}/2)\\}$,
as displayed in Fig. (2).

Figure 2: Some of the 50 grid points $\\{\xi_{k1},\xi_{k2}\\}$ on the surface
of the unit sphere. There are more points on the back, which are invisible
from the present perspective.

The procedure just presented can be generalized to more than just two modes
and will allow us to perform numerical calculations with a finite, presumably
small number of GCS for the BH model. Before we can do so, we have to derive
the relevant equations of motion, however.

##  III Lagrangian formulation of the variational dynamics

Generalizing the variational coherent state Ansatz reviewed in Werther  _et
al._ (2021) by using the GCS discussed above, we are led to the expression

$\displaystyle|\Psi(t)\rangle=\sum_{k=1}^{N}A_{k}(t)|S,\vec{\xi}_{k}(t)\rangle$
(20)

for the wavefunction, where
$|S,\vec{\xi}_{k}(t)\rangle=|S,\xi_{k1}(t),\xi_{k2}(t),\cdots\xi_{kM}(t)\rangle$
and $M$ is the number of modes in the BH Hamiltonian of Eq. (1). The initial
set of basis functions $\\{|S,\vec{\xi}_{k}(0)\rangle\\}$ is the one
constructed in the previous section. Both the expansion coefficients
$\\{A_{k}(t)\\}$ as well as the basis functions are time-dependent. The
multiplicity parameter $N$ determines the number of interacting, complex
parameters whose dynamics is to be determined. It is given by $(M+1)N$ and
will have to be compared to the number of Fock states given in Eq. (11). The
$N$ replicas of condition (10), in principle, reduce the number of independent
real parameters by $N$ but we did not make use of this fact in the numerics to
be presented in the next section (although we monitored that (10) is
fulfilled). We stress that the Ansatz above goes beyond mean-field Ansätze a
la Gutzwiller, used, e.g., in Jaksch  _et al._ (1998).

The time-dependent variational principle can be formulated in terms of the
Lagrangian

$L=\frac{\rm
i}{2}[\langle\Psi|\dot{\Psi}\rangle-\langle\dot{\Psi}|\Psi\rangle]-\langle\Psi(t)|\hat{H}|\Psi(t)\rangle,$
(21)

and is leading to the Euler-Lagrange equations

$\frac{\partial L}{\partial u_{k}^{*}}-\frac{\rm d}{{\rm d}t}\frac{\partial
L}{\partial\dot{u}_{k}^{*}}=0,$ (22)

where $u_{k}$ can be any element of $\\{A_{k},\xi_{k1},\cdots\xi_{kM}\\}$
Kramer and Saraceno (1981). A comparison of the different variational
principles in this context, which are McLachlan, time-dependent and Dirac-
Frenkel is performed in Broekhove  _et al._ (1988) and Lee and Fischer (2014);
Alon  _et al._ (2008) as well as in Werther (2020). From Eq. (20) it is
obvious that the wavefunction Ansatz does not contain the complex conjugates
of any of the parameters and that thus the Cauchy-Riemann conditions are
fulfilled and all three variational principles are equivalent, analogous to
the equivalence worked out in Alon  _et al._ (2008) for MCTDHB.

To proceed, we first need to calculate the expression (suppressing the $S$
dependence of the basis states in a short-hand notation)

$\displaystyle\frac{\rm
i}{2}[\langle\Psi|\dot{\Psi}\rangle-\langle\dot{\Psi}|\Psi\rangle]$
$\displaystyle=\frac{\rm
i}{2}\sum_{k,j=1}^{N}(A_{k}^{*}\dot{A}_{j}-\dot{A}^{*}_{k}A_{j})\langle\vec{\xi}_{k}|\vec{\xi}_{j}\rangle$
$\displaystyle+\frac{\rm
i}{2}S\sum_{k,j=1}^{N}A_{k}^{*}A_{j}\sum_{i=1}^{M}(\xi_{ki}^{*}\dot{\xi}_{ji}-\dot{\xi}_{ki}^{*}\xi_{ji})\langle\vec{\xi^{{}^{\prime}}_{k}}|\vec{\xi^{{}^{\prime}}_{j}}\rangle$
(23)

for the time-derivative. Here we have used the expressions for the right and
left time derivatives from Eqs. (A,A) to arrive at the final result.

Furthermore, from the Hamiltonian for the multimode BH model in Eq. (1), we
get the expectation value

$\displaystyle H$ $\displaystyle=\langle\Psi(t)|\hat{H}|\Psi(t)\rangle$
$\displaystyle=\sum_{k,j=1}^{N}A_{k}^{*}A_{j}\Big{[}-JS\sum_{i=1}^{M-1}(\xi_{ki}^{*}\xi_{j,i+1}+\xi_{k,i+1}^{*}\xi_{ji})\langle\vec{\xi^{{}^{\prime}}_{k}}|\vec{\xi_{j}^{{}^{\prime}}}\rangle+\frac{U}{2}S(S-1)\sum_{i=1}^{M}\xi_{ki}^{*2}\xi_{ji}^{2}\langle\vec{\xi_{k}^{{}^{\prime\prime}}}|\vec{\xi_{j}^{{}^{\prime\prime}}}\rangle$
$\displaystyle+S\frac{K}{2}\sum_{i=1}^{M}(i-j_{0})^{2}\xi_{ki}^{*}\xi_{ji}\langle\vec{\xi^{{}^{\prime}}_{k}}|\vec{\xi^{{}^{\prime}}_{j}}\rangle\Big{]}$
(24)

where we have again used definitions of the $(S-1)$ and $(S-2)$-boson GCS from
Appendix A. Combining Eq. (III) and Eq. (III), we arrive at the Lagrangian for
the multimode BH model

$\displaystyle L$ $\displaystyle=\frac{\rm
i}{2}\sum_{k,j=1}^{N}(A_{k}^{*}\dot{A}_{j}-\dot{A}^{*}_{k}A_{j})\langle\vec{\xi}_{k}|\vec{\xi}_{j}\rangle+\frac{\rm
i}{2}S\sum_{k,j=1}^{N}A_{k}^{*}A_{j}\sum_{i=1}^{M}(\xi_{ki}^{*}\dot{\xi}_{ji}-\dot{\xi}_{ki}^{*}\xi_{ji})\langle\vec{\xi^{{}^{\prime}}_{k}}|\vec{\xi^{{}^{\prime}}_{j}}\rangle$
$\displaystyle-\sum_{k,j=1}^{N}A_{k}^{*}A_{j}\Big{[}-JS\sum_{i=1}^{M-1}(\xi_{ki}^{*}\xi_{j,i+1}+\xi_{k,i+1}^{*}\xi_{ji})\langle\vec{\xi^{{}^{\prime}}_{k}}|\vec{\xi_{j}^{{}^{\prime}}}\rangle+\frac{U}{2}S(S-1)\sum_{i=1}^{M}\xi_{ki}^{*2}\xi_{ji}^{2}\langle\vec{\xi_{k}^{{}^{\prime\prime}}}|\vec{\xi_{j}^{{}^{\prime\prime}}}\rangle$
$\displaystyle+S\frac{K}{2}\sum_{i=1}^{M}(i-j_{0})^{2}\xi_{ki}^{*}\xi_{ji}\langle\vec{\xi^{{}^{\prime}}_{k}}|\vec{\xi^{{}^{\prime}}_{j}}\rangle\Big{]}.$
(25)

This leads to the following derivatives of the Lagrangian with respect to the
coefficients $A_{k}^{*}$ and their time derivatives

$\displaystyle\frac{\partial L}{\partial A^{*}_{k}}$ $\displaystyle=\frac{\rm
i}{2}\sum_{j=1}^{N}\dot{A}_{j}\langle\vec{\xi_{k}}|\vec{\xi_{j}}\rangle+\frac{\rm
i}{2}S\sum_{j=1}^{N}A_{j}\sum_{i=1}^{M}(\xi_{ki}^{*}\dot{\xi}_{ji}-\dot{\xi}_{ki}^{*}\xi_{ji})\langle\vec{\xi^{{}^{\prime}}_{k}}|\vec{\xi^{{}^{\prime}}_{j}}\rangle-\frac{\partial
H}{\partial A_{k}^{*}}$ $\displaystyle\frac{\rm d}{{\rm d}t}\frac{\partial
L}{\partial\dot{A}_{k}^{*}}$ $\displaystyle=-\frac{\rm
i}{2}\sum_{j=1}^{N}\dot{A}_{j}\langle\vec{\xi}_{k}|\vec{\xi}_{j}\rangle-\frac{\rm
i}{2}S\sum_{j=1}^{N}A_{j}\sum_{i=1}^{M}(\xi^{*}_{ki}\dot{\xi}_{ji}+\dot{\xi}_{ki}^{*}\xi_{ji})\langle\vec{\xi_{k}^{{}^{\prime}}}|\vec{\xi^{{}^{\prime}}_{j}}\rangle,$
(26)

where

$\displaystyle\frac{\partial H}{\partial A_{k}^{*}}$
$\displaystyle=\sum_{k,j=1}^{N}A_{j}\Big{[}-JS\sum_{i=1}^{M-1}(\xi_{ki}^{*}\xi_{j,i+1}+\xi_{k,i+1}^{*}\xi_{ji})\langle\vec{\xi^{{}^{\prime}}_{k}}|\vec{\xi_{j}^{{}^{\prime}}}\rangle+\frac{U}{2}S(S-1)\sum_{i=1}^{M}\xi_{ki}^{*2}\xi_{ji}^{2}\langle\vec{\xi_{k}^{{}^{\prime\prime}}}|\vec{\xi_{j}^{{}^{\prime\prime}}}\rangle$
$\displaystyle+S\frac{K}{2}\sum_{i=1}^{M}(i-j_{0})^{2}\xi_{ki}^{*}\xi_{ji}\langle\vec{\xi^{{}^{\prime}}_{k}}|\vec{\xi^{{}^{\prime}}_{j}}\rangle\Big{]}$
(27)

From the Euler Lagrange equation (22), the equation of motion for the
coefficient $A_{j}$

$\displaystyle{\rm
i}\sum_{j=1}^{N}\dot{A}_{j}\langle\vec{\xi_{k}}|\vec{\xi_{j}}\rangle+iS\sum_{j=1}^{N}A_{j}\sum_{i=1}^{M}\xi_{ki}^{*}\dot{\xi}_{ji}\langle\vec{\xi^{{}^{\prime}}_{k}}|\vec{\xi^{{}^{\prime}}_{j}}\rangle-\frac{\partial
H}{\partial A_{k}^{*}}=0$ (28)

can thus be obtained.

Next we switch to the equations of motion for the parameters $\xi_{km}$ of the
SU($M$) coherent states. For the corresponding derivatives of the Lagrangian,
we find

$\displaystyle\frac{\partial L}{\partial\xi^{*}_{km}}$
$\displaystyle=\frac{\rm
i}{2}S\sum_{j=1}^{N}(A_{k}^{*}\dot{A}_{j}-\dot{A}^{*}_{k}A_{j})\xi_{jm}\langle\vec{\xi_{k}^{{}^{\prime}}}|\vec{\xi_{j}^{{}^{\prime}}}\rangle+\frac{\rm
i}{2}S\sum_{j=1}^{N}A_{k}^{*}A_{j}\dot{\xi}_{jm}\langle\vec{\xi_{k}^{{}^{\prime}}}|\vec{\xi_{j}^{{}^{\prime}}}\rangle$
$\displaystyle+\frac{\rm
i}{2}S(S-1)\sum_{j=1}^{N}A_{k}^{*}A_{j}\Big{[}\sum_{i=1}^{M}(\xi_{ki}^{*}\dot{\xi}_{ji}-\dot{\xi}_{ki}^{*}\xi_{ji})\xi_{jm}\langle\vec{\xi_{k}^{{}^{\prime\prime}}}|\vec{\xi_{j}^{{}^{\prime\prime}}}\rangle\Big{]}-\frac{\partial
H}{\partial\xi_{km}^{*}}$ (29) $\displaystyle\frac{\rm d}{{\rm
d}t}\frac{\partial L}{\partial\dot{\xi}_{km}^{*}}$ $\displaystyle=-\frac{\rm
i}{2}S\sum_{j=1}^{N}(\dot{A}_{k}^{*}A_{j}+A_{k}^{*}\dot{A}_{j})\xi_{jm}\langle\vec{\xi_{k}^{{}^{\prime}}}|\vec{\xi_{j}^{{}^{\prime}}}\rangle-\frac{\rm
i}{2}S\sum_{j=1}^{N}A_{k}^{*}A_{j}\dot{\xi}_{jm}\langle\vec{\xi_{k}^{{}^{\prime}}}|\vec{\xi_{j}^{{}^{\prime}}}\rangle$
$\displaystyle-\frac{\rm
i}{2}S(S-1)\sum_{j=1}^{N}A_{k}^{*}A_{j}\sum_{j=1}^{M}(\dot{\xi}_{ki}^{*}\xi_{ji}+\xi_{ki}^{*}\dot{\xi}_{ji})\xi_{jm}\langle\vec{\xi_{k}^{{}^{\prime\prime}}}|\vec{\xi_{j}^{{}^{\prime\prime}}}\rangle$
(30)

with

$\displaystyle\frac{\partial H}{\partial\xi_{km}^{*}}$
$\displaystyle=\sum_{j=1}^{N}A_{k}^{*}A_{j}\Big{[}-JS(\xi_{j,m+1}+\xi_{j,m-1})\langle\vec{\xi_{k}^{{}^{\prime}}}|\vec{\xi_{j}^{{}^{\prime}}}\rangle$
$\displaystyle-
JS(S-1)\sum_{i=1}^{M-1}(\xi_{ki}^{*}\xi_{j,i+1}+\xi_{k,i+1}^{*}\xi_{j,i})\xi_{jm}\langle\vec{\xi_{k}^{{}^{\prime\prime}}}|\vec{\xi_{j}^{{}^{\prime\prime}}}\rangle$
$\displaystyle+US(S-1)\xi_{km}^{*}\xi_{jm}^{2}\langle\vec{\xi_{k}^{{}^{\prime\prime}}}|\vec{\xi_{j}^{{}^{\prime\prime}}}\rangle+\frac{U}{2}S(S-1)(S-2)\sum_{i=1}^{M}\xi_{ki}^{*2}\xi_{ji}^{2}\xi_{jm}\langle\vec{\xi_{k}^{{}^{\prime\prime\prime}}}|\vec{\xi_{j}^{{}^{\prime\prime\prime}}}\rangle$
$\displaystyle+\frac{K}{2}S(m-j_{0})^{2}\xi_{jm}\langle\vec{\xi_{k}^{{}^{\prime}}}|\vec{\xi_{j}^{{}^{\prime}}}\rangle+\frac{K}{2}S(S-1)\sum_{i=1}^{M}(i-j_{0})^{2}\xi_{ki}^{*}\xi_{ji}\xi_{jm}\langle\vec{\xi_{k}^{{}^{\prime\prime}}}|\vec{\xi_{j}^{{}^{\prime\prime}}}\rangle\Big{]}$
(31)

Again using the Euler-Lagrange equation, we finally arrive at the differential
equation

$\displaystyle{\rm
i}S\Bigl{[}\sum_{j=1}^{N}A_{k}^{*}\dot{A}_{j}\xi_{jm}\langle\vec{\xi_{k}^{{}^{\prime}}}|\vec{\xi_{j}^{{}^{\prime}}}\rangle+\sum_{j=1}^{N}A_{k}^{*}A_{j}\dot{\xi}_{jm}\langle\vec{\xi_{k}^{{}^{\prime}}}|\vec{\xi_{j}^{{}^{\prime}}}\rangle+(S-1)\sum_{j=1}^{N}A_{k}^{*}A_{j}\sum_{i=1}^{M}\xi_{ki}^{*}\dot{\xi}_{ji}\xi_{jm}\langle\vec{\xi_{k}^{{}^{\prime\prime}}}|\vec{\xi_{j}^{{}^{\prime\prime}}}\rangle\Bigr{]}$
(32) $\displaystyle-\frac{\partial H}{\partial\xi_{km}^{*}}=0$

for the $\xi_{km}$. Eqs. (28,32) are coupled and nonlinear, as in the more
standard case of Glauber coherent states as basis functions. There are no
exponential function nonlinearities though, because the GCS overlaps are
“simpler” than in the Glauber case, as can be seen in Appendix A.

In Appendix B, we put the nonlinear, coupled, implicit differential equations
for the GCS parameters and the expansion coefficients in matrix form, which
allows us to solve them efficiently numerically.

##  IV Numerical results

As a proof of principle, in the following, we will present numerical results
for the population dynamics of the BH model with $J=1$, based on the complex
grid for SU($M$) coherent states introduced in Sect. II.2 for increasing mode
number $M<S$, leading to an ever increasing size of the Fock state basis. If
possible, numerical results from a calculation in the full Fock space are
compared with the GCS results. As a guide, in Table 1, the different
parameters for which results will be presented are collected, highlighting the
dimension of the Fock space as well as the number of complex parameters,
$(M+1)N$, needed for converged results using the GCS method. We note that
almost converged results can be gained by using much smaller GCS bases than
indicated in the table. Furthermore, for longer time series more basis
functions will be needed and, as shown in Appendix C, closely spaced grids
help to keep the number of basis functions needed for convergence small. For
minimization of the number of basis functions, both, randomization of the
grid, as well as reduction of grid spacing should be performed. This has been
done only for the cases $M=4$ and $M=6$, however.

| 2 | 3 | 4 | 6
---|---|---|---|---
20 |  |  $C^{F}$=231 |  |  $C^{F}$=53130
|  | 200, IV.2 |  | 3500, IV.3
30 |  |  |  $C^{F}$=5456 |
|  |  | 500, App. C, |
50 |  $C^{F}$=51 |  |  |
| 75, IV.1, App. C |  |  |
200 |  $C^{F}$=201 |  |  |
| 243, App. C |  |  |
Table 1: Collection of parameters of numerical studies presented below. In
dependence of mode number $M$ (column index) and particle number $S$ (line
index), the table’s entries give the size of the Fock state basis from Eq.
(11), as well as the number of complex parameters needed for converged GCS
results and where in the text the numerical results can be found. For
$M=2,S=200$, the underlying grid was a diagonal grid (see App. C). For
$M=2,S=50$, the grid in Sec. IV.1 was random with spacing $\sqrt{\pi}$ and the
one in App. C was diagonal but with optimized spacing.

###  IV.1 The two mode case $M=2$

Firstly, for the driven two-mode BH model with a time-dependent tunneling
matrix element studied in Trimborn  _et al._ (2009), the Hamiltonian can be
written as

$\displaystyle
H=-J(t)(a_{1}^{\dagger}a_{2}+a_{2}^{\dagger}a_{1})+\frac{U}{2}(a_{1}^{{\dagger}2}a_{1}^{2}+a_{2}^{{\dagger}2}a_{2}^{2})+\frac{K}{2}\sum_{j=1}(j-j_{0})^{2}a_{j}^{\dagger}a_{j}$
(33)

where $J(t)=J_{0}+J_{1}\cos(\omega t)$. The initial state

$|\Psi(0)\rangle=|S=50,\xi_{1}=-\sqrt{0.7},\xi_{2}=\sqrt{0.3}\rangle$ (34)

used for the propagation is the same as in Trimborn  _et al._ (2009).

By means of the SU($2$) coherent states originating from the complex grid of
Glauber coherent with distance $\sqrt{\pi}$, we employ the Ansatz

$|\Psi(t)\rangle=\sum_{k=1}^{N}A_{k}(t)|S=50,\xi_{k1}(t),\xi_{k2}(t)\rangle$
(35)

for the wavefunction that contains 3$N$ complex parameters. The random grid
one is constructed such that it contains the initial state in the “middle”, as
indicated in Fig. 1.

Converged results for the population (normalized to the total number of
particles) of the first mode as a function of time were obtained by using
$N=$25 basis states and even 15 states lead to almost perfect results as is
shown in Fig. (3). These results agree with the ones shown in Fig. 12(c) in
Trimborn  _et al._ (2009). We stress that for trivial multiplicity $N=1$, the
result is only reasonable for a very short time period and thereafter looses
its predictive power, as can be seen in Fig. (3).

There is no big difference in the number of complex parameters that have been
used for the converged results, compared to the size of the Fock space basis,
which here is just 51, however. Therefore let us investigate increasing mode
number in the following.

Figure 3: Dynamics of the population of the first mode of a two mode BH model
with different numbers of SU(2) basis states. The particle number is $S=50$,
the hopping parameter is $J_{0}=1$, the on-site interaction energy is
$U=0.1J_{0}$ and $K=0$. The driving frequency and strength are
$\omega=2\pi/J_{0}$ and $J_{1}=0.5J_{0}$. Different lines are results for
different numbers of GCS: solid: $N=1$, dashed: $N=15$, dotted: $N=25$, dash-
dotted: $N=50$. Stars denote Fock space results.

###  IV.2 The three mode case $M=3$

For three modes occupied by 20 bosons, we first choose the superposition of
Fock states

$|\Psi(0)\rangle=|S=20,\xi_{1}=\frac{\sqrt{2}}{2},\xi_{2}=i\frac{\sqrt{2}}{2},\xi_{3}=0\rangle$
(36)

as the initial state. In this case, our numerical results in Fig. (4) show
that convergence is achieved with 50 basis functions, i.e. 200 complex
parameters, a little less than the 231 complex-valued amplitudes one would
need for an expansion in Fock states in the present case.

Figure 4: Dynamics of the population of the first mode for three mode case for
initial state (36). The on-site interaction energy is $U=0.1J$, whereas $K=0$.
Results for different number of basis functions (solid: $N=50$, dashed
$N=100$, dash-dotted $N=200$) are almost indistinguishable. Stars denote Fock
space results.

The parameter combination $\Lambda={US/J}$ is frequently applied to
distinguish the dynamical features of the BH model. For $\Lambda<1$ the
dynamics is located in Rabi regime, while $1<\Lambda\ll S^{2}$ and $\Lambda\gg
S^{2}$ represent the so-called Josephson and Fock regime, respectively Leggett
(2001); Simon and Strunz (2014). In Figs. 5 and 6, we show the accuracy of the
variational dynamics for the initial Fock state

$|\Psi(0)\rangle=|S=20,1,0,0\rangle$ (37)

by comparing with the exact numerical results (gained by an expansion of the
wavefuncion in terms of Fock states) for 3 modes and 20 bosons. It turns out
that the 50 basis functions that were found sufficient in Fig. 4 are still
sufficient in the Rabi and Josephson cases. The pure Fock regime case was too
demanding numerically and is therefore not considered. In the Rabi case
displayed in Fig. 5 almost undamped sinusoidal oscillations are observed. The
oscillations become more complex and damped in the Josephson case shown in
Fig. 6. In this case, also the accuracy of ode45 MATLAB (2020) used for the
solution of the differential equations had to be promoted beyond the default
value.

Figure 5: Dynamics of the population of all three modes. The initial state is
the Fock state $|20,1,0,0\rangle$. $U=0.03J$ and $\Lambda=0.6$ is in the Rabi
regime. The exact results are displayed by the different markers (crosses:
$\langle N_{1}\rangle/S$, triangles: $\langle N_{2}\rangle/S$, circles:
$\langle N_{3}\rangle/S$), while the corresponding lines without markers are
calculated by variational dynamics using 50 basis functions whose distance is
$\sqrt{\pi}$.
Figure 6: Dynamics of the population of all three modes. The initial state is
the Fock state $|20,1,0,0\rangle$. $U=0.2J$ and $\Lambda=4$ is in the
Josephson regime. The exact results are displayed by the different markers
(crosses: $\langle N_{1}\rangle/S$, triangles: $\langle N_{2}\rangle/S$,
circles: $\langle N_{3}\rangle/S$), while the corresponding lines without
markers are calculated by the variational dynamics using 50 basis functions
whose distance is $\sqrt{\pi}$.

###  IV.3 Beyond $M=3$

To highlight the power of the proposed method we now investigate the more
complicated situation where $20$ bosons occupy $6$ modes. Then the Hilbert
space consists of ${25\choose 5}=53130$ Fock states and with a standard
desktop computer, it will be difficult to propagate in this large Hilbert
space using the full Fock state basis.

If we employ the SU$(M)$ states originating from the complex grid with the von
Neumann spacing $\sqrt{\pi}$, as presented in Fig. (7), we find that only
several hundreds of GCS are needed for converged results at longer times. This
means we need $(M+1)\times 800=5600$ complex valued parameters in our Ansatz
(20), in order to converge the results for all times shown. Especially at the
extrema of the curve, a large number of basis functions is needed.

Figure 7: The dynamics of population of the first mode for the case of six
modes. The initial state is
$|S=20,\xi_{1}=\frac{\sqrt{2}}{2},\xi_{2}=i\frac{\sqrt{2}}{2},\xi_{3}=\cdots\xi_{6}=0\rangle$.
The on-site interaction energy is $U=0.1J$ and $K=0$. Solid line: 600 basis
functions, dashed line: 700 basis functions, dashed-dotted line: 800 basis
functions, dotted line: 1000 basis functions. The spacing of the underlying
Glauber CS grid was $\sqrt{\pi}$.

However, if we decrease the distance of the complex grid, the results can be
converged faster. In Fig. (8), we show that the underlying grid with the
smaller distance of $\sqrt{\pi}/32$ between the grid points allows for
convergence already with 500 basis functions, whereas the “standard” grid
needed 800 states to converge the results at specific later times. An in-depth
study of the usage of denser underlying grids is given in Appendix C.

Figure 8: Convergence of the results for a distance of $\frac{\sqrt{\pi}}{32}$
of the underlying Glauber CS grid. Other parameters (apart from the number of
basis functions) as in Fig. (7)

These numerical results show that a GCS with only $M$ complex parameters
according to Eq. (9) and a corresponding variational Ansatz (20) with $(M+1)N$
parameters allows for a faithful description of the dynamics under the BH
Hamiltonian. The number of independent complex parameters that we had to
employ is more than one order of magnitude smaller than the dimension of the
complete Fock space.

##  V Conclusions and Outlook

We have shown that a numerically favorable treatment of BH dynamics in 1D is
possible by using a variational approach based on GCS. In contrast to
established mean-field versions of the theory, here we have employed a
numerically complete set of SU($M$) GCS as basis states in the expansion of
the initial wavefunction. We have built an overcomplete set of basis function
by using the expansion of the SU($M$) states in terms of Glauber coherent
states for whom it is well-established how a numerically complete set of basis
states has to be chosen.

The equations of motion of the GCS parameters as well as of the expansion
coefficients have been derived from the time-dependent variational principle
in its Lagrangian form. The central idea is that the number of GCS parameters
times the number of basis functions is much smaller than the number of basis
functions needed in an unbiased expansion of the wavefunction in terms of Fock
states. In this respect the method is closely related to the multi-
configuration time-dependent Hartree method for bosons Alon  _et al._ (2008)
as well as the variational approach in Rau  _et al._ (2010a, b), which are
both, however, not based on GCS.

That the method indeed works for considerably smaller numbers of basis
functions than in the case of an expansion in terms of Fock states has shown
to be true for large particle number $S$ and increasing mode number $M<S$ for
several examples in different parameter ranges. Especially for larger mode
number, the use of small spacing (considerably smaller than the von Neumann
spacing of $\sqrt{\pi}$) in the underlying Glauber CS grid turned out to be
necessary for the use of a relatively small number of GCS basis functions. In
this way we could report converged results using more than one order of
magnitude less parameters than would have been necessary for an expansion in
terms of Fock states. We stress that we did not encounter any problems that
had to be solved with the recently introduced apoptosis procedure Werther and
Grossmann (2020a). Regularization also used in Werther and Grossmann (2020a)
and well known from other approaches like MCTDH Manthe  _et al._ (1992) was
necessary, however.

In the future, we plan to investigate also larger mode numbers and small
particle numbers that are relevant in thermalization studies of cold atoms in
optical superlattices Trotzky  _et al._ (2012), as well as other particle
number conserving bosonic lattice systems with the proposed method. If its
fortunate scaling properties persist, new territory in parameter space may
become explorable and we are especially interested in the case where $U$ is
much larger than $J$. Furthermore, the application of the proposed methodology
to BH models in more than 1D is planned. It is well known that the
entanglement entropy of the ground state in lattice models fulfills so-called
area laws Eisert  _et al._ (2010). The entanglement growth in the course of
time evolution poses a serious problem for other numerical propagation
techniques like matrix product state based methods, however Schollwöck (2011).
It remains to be explored how well the method proposed herein can cope with
these problems in higher dimensions.

Acknowledgments: The authors would like to thank Valentin Link for
enlightening discussions. FG would like to thank the Deutsche
Forschungsgemeinschaft for financial support under grant GR 1210/8-1.

##  Appendix A Some properties of SU($M$) coherent states

In this appendix, we review some computationally helpful formulae along the
lines of Buonsante and Penna (2008), which are needed in the main text.
Firstly, the commutation between the annihilation operator and the collective
creation operator which is applied to generate an SU($M$) coherent state is
given by

$\displaystyle[\hat{a}_{i},\Big{(}\sum_{j=1}^{M}\xi_{j}\hat{a}_{j}^{\dagger}\Big{)}^{S}]=S\xi_{i}\Big{(}\sum_{j=1}^{M}\xi_{j}\hat{a}_{j}^{\dagger}\Big{)}^{S-1}.$
(38)

Secondly, by defining two collective operators

$\displaystyle\hat{A}^{\dagger}=\sum_{i=1}^{M}\xi_{i}\hat{a}_{i}^{\dagger},\indent\hat{B}^{\dagger}=\sum_{i=1}^{M}\eta_{i}\hat{a}_{i}^{\dagger},$
(39)

two different SU($M$) coherent states can be generated via

$\displaystyle|\vec{\xi}\,\rangle=\frac{1}{\sqrt{S!}}(\hat{A}^{\dagger})^{S}|0\rangle,\indent|\vec{\eta}\,\rangle=\frac{1}{\sqrt{S!}}(\hat{B}^{\dagger})^{S}|0\rangle.$
(40)

The vectorized parameter contains the parameters of all modes (here we do not
consider doubly indexed parameters). The inner product of the above two states
is then given by

$\displaystyle\langle\vec{\eta}|\vec{\xi}\,\rangle$
$\displaystyle=\frac{1}{S!}\langle
0|\Big{(}\sum_{i=1}^{M}\eta^{*}_{i}\hat{a}_{i}\Big{)}^{S}\Big{(}\sum_{j=1}^{M}\xi_{j}\hat{a}_{j}^{\dagger}\Big{)}^{S}|0\rangle$
$\displaystyle=\frac{1}{S!}\langle
0|\sum_{n_{1}+n_{2}+\cdots=S}\frac{S!}{n_{1}!n_{2}!\cdots}[(\eta_{1}^{*}\hat{a}_{1})^{n_{1}}(\eta_{2}^{*}\hat{a}_{2})^{n_{2}}\cdots]$
$\displaystyle\sum_{m_{1}+m_{2}+\cdots=S}\frac{S!}{m_{1}!m_{2}!\cdots}[(\xi_{1}\hat{a}_{1}^{\dagger})^{m_{1}}(\xi_{2}\hat{a}_{2}^{\dagger})^{m_{2}}\cdots]|0\rangle$
$\displaystyle=\frac{1}{S!}\Big{(}\sum_{n_{1}+n_{2}+\cdots=S}\frac{S!}{\sqrt{n_{1}!n_{2}!\cdots}}\eta_{1}^{*n_{1}}\eta_{2}^{*n_{2}}\cdots\langle\vec{n}|\Big{)}$
$\displaystyle\sum_{m_{1}+m_{2}+\cdots=S}\frac{S!}{\sqrt{m_{1}!m_{2}!\cdots}}\xi_{1}^{m_{1}}\xi_{2}^{m_{2}}\cdots|\vec{m}\rangle$
$\displaystyle=\sum_{m_{1}+m_{2}+\cdots=S}\frac{S!}{m_{1}!m_{2}!\cdots}(\eta^{*}_{1}\xi_{1})^{m_{1}}(\eta_{2}^{*}\xi_{2})^{m_{2}}\cdots$
$\displaystyle=\Big{(}\sum_{i=1}^{M}\eta_{i}^{*}\xi_{i}\Big{)}^{S},$ (41)

where we have used the general binomial theorem

$\displaystyle(x_{1}+x_{2}+\cdots+x_{n})^{k}=\sum_{a_{1}+a_{2}+\cdots+a_{n}=k}\frac{k!}{a_{1}!a_{2}!\cdots
a_{n}!}x_{1}^{a_{1}}x_{2}^{a_{2}}\cdots x_{n}^{a_{n}}.$ (42)

The above result (41) is quite different from the corresponding property for
Glauber coherent states, which involves exponential functions.

Using Eq. (38), we can now calculate the action of the annihilation operator
on the SU($M$) coherent state via

$\displaystyle\hat{a}_{i}|\vec{\xi}\,\rangle$
$\displaystyle=\hat{a}_{i}\frac{1}{\sqrt{S!}}\Big{(}\sum_{j=1}^{M}\xi_{j}\hat{a}_{j}^{\dagger}\Big{)}^{S}|0\rangle$
$\displaystyle=\frac{1}{\sqrt{S!}}[\Big{(}\sum_{j=1}^{M}\xi_{j}\hat{a}_{j}^{\dagger}\Big{)}^{S}a_{i}+S\xi_{i}\Big{(}\sum_{j=1}^{M}\xi_{j}\hat{a}_{j}^{\dagger}\Big{)}^{S-1}]|0\rangle$
$\displaystyle=\sqrt{S}\xi_{i}|\vec{\xi^{{}^{\prime}}}\rangle$ (43)

where we have used the action of the annihilation operator on the ground
state, see Eq. (2) as well as the definition

$|\vec{\xi^{{}^{\prime}}}\rangle=\frac{1}{\sqrt{(S-1)!}}\Big{(}\sum_{j=1}^{M}\xi_{j}\hat{a}_{j}^{\dagger}\Big{)}^{S-1}|0\rangle$
(44)

of the $(S-1)$-boson GCS.

Furthermore, for $|\eta\rangle$ and $|\xi\rangle$, from Eq. (A) we get

$\displaystyle\langle\vec{\eta}|\hat{a}_{j}^{\dagger}\hat{a}_{k}|\vec{\xi}\,\rangle=S\eta^{*}_{j}\xi_{k}\langle\vec{\eta^{{}^{\prime}}}|\vec{\xi^{{}^{\prime}}}\rangle$
(45)

where the inner product
$\langle\vec{\eta^{{}^{\prime}}}|\vec{\xi^{{}^{\prime}}}\rangle$ is

$\displaystyle\langle\vec{\eta^{{}^{\prime}}}|\vec{\xi^{{}^{\prime}}}\rangle=\Big{(}\sum_{i=1}^{M}\eta_{i}^{*}\xi_{i}\Big{)}^{S-1}.$
(46)

In the main text as well as below, we will also refer to the following two
inner products

$\displaystyle\langle\vec{\eta^{{}^{\prime\prime}}}|\vec{\xi^{{}^{\prime\prime}}}\rangle$
$\displaystyle=\Big{(}\sum_{i=1}^{M}\eta_{i}^{*}\xi_{i}\Big{)}^{S-2},$ (47)
$\displaystyle\langle\vec{\eta^{{}^{\prime\prime\prime}}}|\vec{\xi^{{}^{\prime\prime\prime}}}\rangle$
$\displaystyle=\Big{(}\sum_{i=1}^{M}\eta_{i}^{*}\xi_{i}\Big{)}^{S-3},$ (48)

of the $(S-2)$ and the $(S-3)$ GCS.

Using the chain rule, the matrix element of the right time-derivative now is

$\displaystyle\langle\vec{\eta}|\overrightarrow{\partial_{t}}|\vec{\xi}\,\rangle$
$\displaystyle=\langle\vec{\eta}|\sum_{i=1}^{M}\dot{\xi_{i}}\partial_{\xi_{i}}|\vec{\xi}\,\rangle$
$\displaystyle=\frac{S}{\sqrt{S!}}\langle\vec{\eta}|\sum_{i=1}^{M}(\dot{\xi}_{i}\hat{a}_{i}^{\dagger})(\hat{A}^{\dagger})^{S-1}|0\rangle$
$\displaystyle=\frac{S}{S!}\langle
0|\hat{B}^{S-1}\hat{B}\sum_{i=1}^{M}(\dot{\xi}_{i}\hat{a}_{i}^{\dagger})(\hat{A}^{\dagger})^{S-1}|0\rangle$
$\displaystyle=\langle\vec{\eta^{{}^{\prime}}}|\sum_{i=1}^{M}(\dot{\xi}_{i}\hat{B}\hat{a}_{i}^{\dagger})|\vec{\xi^{{}^{\prime}}}\rangle$
$\displaystyle=\langle\vec{\eta^{{}^{\prime}}}|\sum_{i=1}^{M}\dot{\xi}_{i}(\eta_{i}^{*}+\hat{a}_{i}^{\dagger}\hat{B})|\vec{\xi^{{}^{\prime}}}\rangle$
$\displaystyle=\sum_{i=1}^{M}(\dot{\xi}_{i}\eta^{*}_{i})\langle\vec{\eta^{{}^{\prime}}}|\vec{\xi^{{}^{\prime}}}\rangle+\sum_{i,j=1}^{M}\dot{\xi}_{i}\eta^{*}_{j}\langle\vec{\eta^{{}^{\prime}}}|\hat{a}_{i}^{\dagger}\hat{a}_{j}|\vec{\xi^{{}^{\prime}}}\rangle$
$\displaystyle=\sum_{i=1}^{M}(\dot{\xi}_{i}\eta^{*}_{i})\langle\vec{\eta^{{}^{\prime}}}|\vec{\xi^{{}^{\prime}}}\rangle+(S-1)\sum_{i,j=1}^{M}\dot{\xi}_{i}\eta^{*}_{j}\eta_{i}^{*}\xi_{j}\langle\vec{\eta^{{}^{\prime\prime}}}|\vec{\xi^{{}^{\prime\prime}}}\rangle$
$\displaystyle=S\sum_{i=1}^{M}(\dot{\xi}_{i}\eta^{*}_{i})\langle\vec{\eta^{{}^{\prime}}}|\vec{\xi^{{}^{\prime}}}\rangle$
(49)

In the fifth line of the above equation, we have used the result of Eq. (38)
for $S=1$, and the relation
$\sum_{j=1}^{M}\xi_{j}\eta_{j}^{*}\langle\vec{\eta^{{}^{\prime\prime}}}|\vec{\xi^{{}^{\prime\prime}}}\rangle=\langle\vec{\eta^{{}^{\prime}}}|\vec{\xi^{{}^{\prime}}}\rangle$
which follows from Eqs. (46,47) is also used to get the result of the last
line.

Similarly, for the left time-derivative we find

$\displaystyle\langle\vec{\eta}|\overleftarrow{\partial_{t}}|\vec{\xi}\,\rangle$
$\displaystyle=\langle\vec{\eta}|\sum_{i=1}^{M}\dot{\eta^{*}_{i}}\partial_{\eta^{*}_{i}}|\vec{\xi}\rangle$
$\displaystyle=\langle\vec{\eta^{{}^{\prime}}}|\sum_{i=1}^{M}(\dot{\eta_{i}^{*}}\hat{a}_{i}\hat{A}^{\dagger})|\vec{\xi^{{}^{\prime}}}\rangle$
$\displaystyle=\sum_{i=1}^{M}(\dot{\eta^{*}_{i}}\xi_{i})\langle\vec{\eta^{{}^{\prime}}}|\vec{\xi^{{}^{\prime}}}\rangle+(S-1)\sum_{i,j=1}^{M}(\dot{\eta^{*}_{i}}\xi_{j}\eta^{*}_{j}\xi_{i})\langle\vec{\eta^{{}^{\prime\prime}}}|\vec{\xi^{{}^{\prime\prime}}}\rangle$
$\displaystyle=S\sum_{i=1}^{M}(\dot{\eta^{*}_{i}}\xi_{i})\langle\vec{\eta^{{}^{\prime}}}|\vec{\xi^{{}^{\prime}}}\rangle$
(50)

Both results will be used in the main text.

##  Appendix B Matrix form of the variational equations

Combining the equation Eq. (28) and Eq. (32), we get a compact and scalable
matrix equation for the vector containing all coefficients $\\{A_{k}\\}$ and
GCS parameters $\\{\xi_{km}\\}$

$\left(\begin{array}[]{cc}\bm{X}&\bm{Y}\\\
\bm{Y}^{\dagger}&\bm{Z}\end{array}\right)\left(\begin{array}[]{cc}\bm{\dot{A}}\\\
\bm{\dot{\xi}}\end{array}\right)=-{\rm
i}\left(\begin{array}[]{cc}\bm{R}_{1}\\\ \bm{R}_{2}\end{array}\right),$ (51)

analogous to the procedure lined out in the appendix of Werther and Grossmann
(2020a) for Glauber coherent states.

The block matrices are given by

$\displaystyle\bm{X}_{kj}$
$\displaystyle=\langle\vec{\xi_{k}}|\vec{\xi_{j}}\rangle,$ (52)
$\displaystyle\bm{Y}$
$\displaystyle=S(\bm{\xi}_{1}^{*},\bm{\xi}_{2}^{*},\cdots,\bm{\xi}_{M}^{*})\otimes\bm{A}^{T}\circ(\bm{1}_{1\times
M}\otimes\bm{X}^{{}^{\prime}}),$ (53) $\displaystyle\bm{Z}$
$\displaystyle=\bm{1}_{M\times M}\otimes\bm{\rho}\circ\bm{F},$ (54)

where the vector $\vec{\xi_{k}}$ is now indexed by the basis function
discretization index and the vector $\bm{\xi}_{m}$, to be defined below, is
indexed by the mode index. Furthermore,

$\displaystyle\bm{F}_{ij}$
$\displaystyle=S(S-1)\bm{X^{{}^{\prime\prime}}}\circ({\bm{\xi}}_{j}^{*}\cdot{\bm{\xi}}_{i}^{T})(i\neq
j),$ (55) $\displaystyle\bm{F}_{ii}$
$\displaystyle=S\bm{X^{{}^{\prime}}}+S(S-1)\bm{X^{{}^{\prime\prime}}}\circ({\bm{\xi}}_{i}^{*}\cdot{\bm{\xi}}_{i}^{T}),$
(56)

and where $\bm{1}_{m\times n}$ is an $m\times n$ matrix which only consists of
ones, and
$\bm{X}^{{}^{\prime}}_{kj}=\langle\vec{\xi^{\prime}_{k}}|\vec{\xi^{\prime}_{j}}\rangle$
and
$\bm{X}^{{}^{\prime\prime}}_{kj}=\langle\vec{\xi^{\prime\prime}_{k}}|\vec{\xi^{\prime\prime}_{j}}\rangle$
are overlaps of $(S-1)$ and $(S-2)$-boson GCS from the previous appendix,
respectively, whereas $\bm{\rho}_{kj}=A_{k}^{*}A_{j}$. Furthermore, $\otimes$
denotes the tensor-product, whereas $\circ$ denoted the Hadamard-product
(element-wise multiplication) and $\cdot$ denotes the standard scalar product.

Furthermore, the vectors are defined as

$\bm{\dot{A}}=\left(\begin{array}[]{ccc}\dot{A}_{1}\\\ \dot{A}_{2}\\\
\vdots\\\
\dot{A}_{N}\end{array}\right),\bm{\dot{\xi}}=\left(\begin{array}[]{ccc}\bm{\dot{\xi}_{1}}\\\
\bm{\dot{\xi}_{2}}\\\ \vdots\\\
\bm{\dot{\xi}_{M}}\end{array}\right),\bm{\xi_{m}}=\left(\begin{array}[]{cccc}\xi_{1m}\\\
\xi_{2m}\\\ \vdots\\\
\xi_{Nm}\end{array}\right),\bm{R_{1}}=\left(\begin{array}[]{ccc}\frac{\partial
H}{\partial A_{1}^{*}}\\\ \frac{\partial H}{\partial A_{2}^{*}}\\\ \vdots\\\
\frac{\partial H}{\partial
A_{N}^{*}}\end{array}\right),\bm{R_{2}}=\left(\begin{array}[]{cccc}\frac{\partial
H}{\partial\xi_{11}^{*}}\\\ \frac{\partial H}{\partial\xi_{21}^{*}}\\\
\vdots\\\ \frac{\partial H}{\partial\xi_{1M}^{*}}\\\ \frac{\partial
H}{\partial\xi_{2M}^{*}}\\\ \vdots\end{array}\right)$ (57)

with $H=\langle\Psi|H|\Psi\rangle$, the expectation of the Hamiltonian given
in Eq. (III).

The vectors on the right hand side of equation Eq. (51) are

$\displaystyle\bm{R}_{1}$
$\displaystyle=-JS\sum_{i=1}^{M-1}\big{[}\bm{X}^{{}^{\prime}}\cdot(\bm{A}\circ\bm{\xi_{i+1}})\circ\bm{\xi}_{i}^{*}+\bm{X}^{{}^{\prime}}\cdot(\bm{A}\circ\bm{\xi}_{i})\circ\bm{\xi}_{i+1}^{*}\big{]}$
$\displaystyle+\frac{U}{2}S(S-1)\sum_{i=1}^{M}\bm{X}^{{}^{\prime\prime}}\cdot(\bm{A}\circ\bm{\xi}_{i}\circ\bm{\xi}_{i})\circ\bm{\xi}_{i}^{*}\circ\bm{\xi}_{i}^{*}+\frac{K}{2}S\sum_{i=1}^{M}(i-j_{0})^{2}\bm{X}^{{}^{\prime}}\cdot(\bm{A}\circ\bm{\xi}_{i})\circ\bm{\xi}_{i}^{*}$
(58)

as well as

$\displaystyle\frac{\partial H}{\partial\bm{\xi}_{m}^{*}}$
$\displaystyle=-JS\bm{\rho}\circ\bm{X}^{{}^{\prime}}\cdot(\bm{\xi}_{m+1}+\bm{\xi}_{m-1})$
$\displaystyle-
JS(S-1)\sum_{i=1}^{M-1}\big{[}\bm{\rho}\circ\bm{X}^{{}^{\prime\prime}}\cdot(\bm{\xi}_{m}\circ\bm{\xi}_{i})\circ\bm{\xi}^{*}_{i+1}+\bm{\rho}\circ\bm{X}^{{}^{\prime\prime}}\cdot(\bm{\xi}_{m}\circ\bm{\xi}_{i+1})\circ\bm{\xi}_{i}^{*}\big{]}$
$\displaystyle+US(S-1)\bm{\rho}\circ\bm{X}^{{}^{\prime\prime}}\cdot(\bm{\xi}_{m}\circ\bm{\xi}_{m})\circ\bm{\xi}_{m}^{*}$
$\displaystyle+\frac{U}{2}S(S-1)(S-2)\sum_{i=1}^{M}\big{[}\bm{\rho}\circ\bm{X}^{{}^{\prime\prime\prime}}\cdot(\bm{\xi}_{i}\circ\bm{\xi}_{i}\circ\bm{\xi}_{m})\circ(\bm{\xi}_{i}^{*}\circ\bm{\xi}_{i}^{*})\big{]}$
$\displaystyle+\frac{K}{2}S(m-j_{0})^{2}\bm{\rho}\circ\bm{X}^{{}^{\prime}}\cdot\bm{\xi}_{m}+\frac{K}{2}S(S-1)\sum_{i=1}^{M}\left[(i-j_{0})^{2}\bm{\rho}\circ\bm{X}^{{}^{\prime\prime}}\cdot(\bm{\xi}_{m}\circ\bm{\xi}_{i})\circ\bm{\xi}_{i}^{*}\right]$
(59)

for use in the vector ${\bm{R}}_{2}$. Please keep in mind that $m+1\leq
M,m-1\geq 1$.

##  Appendix C Detailed convergence study with respect to grid spacing

To prove the advantage of small spacings of complex grids for getting
converged results, for most of this appendix we choose the diagonal of the
combination of the rectangular grids rather than the random choice of points
from the rectangular grids which have been used in the main text. In diagonal
grids, the parameters of the SU($M$) coherent state lie in the exact same
positions for every complex subgrid, for example,

$\\{\xi_{mn,1},\xi_{mn,2}\\}=\\{\frac{z_{mn,1}}{\sqrt{|z_{mn,1}|^{2}+|z_{mn,2}|^{2}}},\frac{z_{mn,2}}{\sqrt{|z_{mn,1}|^{2}+|z_{mn,2}|^{2}}}\\},$
(60)

in the case of 2 modes, where $m$, $n$ are the same(!) complex grid indices.
We remind of the fact that in the random choice of the two-mode grid, a random
combination of the points on the grid of the first mode with the points of the
grid of the second mode is allowed! Although using the diagonal grids is less
optimal than the random ones (as will become obvious below), according to our
numerical results, it can help us avoid the influence of randomness and
compare the effects of different spacings directly at the same level. In the
following, we will show some results for normalized population dynamics for
different initial states with different photon and mode numbers.

Figure 9: The dynamics of the population of the first mode in the two-mode
case with 50 photons. The driving frequency and strength are
$\omega=2\pi/J_{0}$ and $J_{1}=0.5J_{0}$, and the on-site interaction energy
$U=0.1J$. The number of SU($M$) coherent states is 25. Different lines are
results for different spacings of the complex grid. Solid: $\sqrt{\pi}$,
dotted: $\sqrt{\pi}/3$, dash: $\sqrt{\pi}/4$, dash-dotted: $\sqrt{\pi}/8$.

For the two-mode case, firstly, as in the main text, the initial state is
chosen as $|50,-\sqrt{0.7},\sqrt{0.3}\rangle$ and the parameters of the
Hamiltonian are also consistent with the ones used in Sect. IV.1. We
constructed two identical complex 5$\times$5 grids, and the number of SU($M$)
we used thus is just $N=25$ (and not 625, as it would be if we would allow all
combinations of the rectangular grid points). Fig. (9) reveals that small
spacings allow to arrive at the converged result more quickly. In the present
case, $\sqrt{\pi}/4$ is enough to reproduce the exact result. Please note that
the result with 25 basis functions and the largest spacing of $\sqrt{\pi}$ is
much worse than the one presented in Fig. 3, where a random grid was used. For
the second, more demanding initial state $|200,-\sqrt{0.7},\sqrt{0.3}\rangle$
with 200 photons, the number of grid points is increased to 81. Fig. (10)
shows that the spacing with $\sqrt{\pi}/4$ does not yet work well and that for
this initial condition $\sqrt{\pi}/8$ is a better choice.

Figure 10: The dynamics of the population of the first mode in the two-mode
case with 200 photons. The driving frequency and strength are
$\omega=2\pi/J_{0}$ and $J_{1}=0.5J_{0}$, and the on-site interaction energy
$U=0.1J$. The number of SU($M$) coherent states is 81. Different lines are
results for different spacings of the complex grid: solid: $\sqrt{\pi}$,
dotted: $\sqrt{\pi}/4$, dash: $\sqrt{\pi}/8$, dash-dotted: $\sqrt{\pi}/16$.

Next we extend our results to the four-mode case with a constant in time
hopping parameter $J=1$. In Fig. (11) results are shown for the initial state
$|30,-\sqrt{0.7},\sqrt{0.3},0,0\rangle$. Although it has less photons than the
in two-mode case above, the spacing$\sqrt{\pi}/4$ which performs well in Fig.
(9) does not give rise to converged results with 169 basis functions, and this
kind of deviation for $\sqrt{\pi}/4$ also occurs for other states with larger
numbers of photons or other initial state parameters (not shown). Thus the
number of modes seems to be decisive, when choosing the optimal grid spacing.

In Fig. (12), we present a study of the number of basis functions which is
needed for getting converged results by optimizing the underlying grid
spacing. This time, however, in order for faster convergence, we go back to
the case of random grids employed in the main text. All results initially
coincide to within line thickness and only for longer times the small spacing
turns out to be advantageous. Although small spacings make the calculation
generally more efficient, decreasing the spacing indefinitely is not an option
for improvement and we found an optimum value of $\sqrt{\pi}/32$ in the
present case. The even smaller spacing $\sqrt{\pi}/64$ will not lead to a
further promotion compared with $\sqrt{\pi}/32$ (the results are even slightly
worse). Our final conclusion is that large numbers of photons and modes will
increase the dimension of Hilbert space dramatically (see Eq. (11)) and make
the dynamical process more complicated. As we have shown here, constructing
complex grids with optimal spacings (smaller than $\sqrt{\pi}$) can help
speeding up the convergence of numerical calculations, and for more
complicated systems smaller spacing is needed for feasibility of the
calculation, but there is no gain in decreasing the spacing indefinitely.

Figure 11: The dynamics of the population of the first mode in the four-mode
case with 30 photons. The on-site interaction energy $U=0.1J$. The number of
SU($M$) coherent states used is 169. Different lines are results for different
spacings of the complex grid: solid: $\sqrt{\pi}$, dotted: $\sqrt{\pi}/4$,
dash: $\sqrt{\pi}/16$, dash-dotted: $\sqrt{\pi}/32$. The crosses display the
exact results using the full Fock state basis.
Figure 12: Convergence of the population dynamics in the four-mode case with
30 photons for different spacings and basis size $N$. The on-site interaction
energy $U=0.1J$. The different curves (solid: $\sqrt{\pi}$, dotted:
$\sqrt{\pi}/4$, dash: $\sqrt{\pi}/16$, dash-dotted: $\sqrt{\pi}/32$, star:
$\sqrt{\pi}/64$, circles: exact Fock results) coincide to within line
thickness for most of the time interval shown.

## References

  * Brif and Mann (1998) C. Brif and A. Mann, Journal of Physics A: Mathematical and General 31, L9 (1998).
  * Brif and Mann (1999) C. Brif and A. Mann, Phys. Rev. A 59, 971 (1999).
  * Trimborn  _et al._ (2008) F. Trimborn, D. Witthaut, and H. J. Korsch, Phys. Rev. A 77, 043631 (2008).
  * Trimborn  _et al._ (2009) F. Trimborn, D. Witthaut, and H. J. Korsch, Phys. Rev. A 79, 013608 (2009).
  * Zhang  _et al._ (1990) W.-M. Zhang, D. H. Feng, and R. Gilmore, Rev. Mod. Phys. 62, 867 (1990).
  * Gardiner and Zoller (2004) C. W. Gardiner and P. Zoller,  _Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics_ , 3rd ed. (Springer-Verlag, Berlin, 2004).
  * Perelomov (1986) A. Perelomov,  _Generalized Coherent States and Their Applications_ (Springer-Verlag, Berlin, 1986).
  * Buonsante and Penna (2008) P. Buonsante and V. Penna, Journal of Physics A: Mathematical and Theoretical 41, 175301 (2008).
  * Arecchi  _et al._ (1972) F. T. Arecchi, E. Courtens, R. Gilmore, and H. Thomas, Phys. Rev. A 6, 2211 (1972).
  * Klimov  _et al._ (2017) A. B. Klimov, J. L. Romero, and H. de Guise, Journal of Physics A: Mathematical and Theoretical 50, 323001 (2017).
  * Merkel  _et al._ (2020) K. Merkel, V. Link, K. Luoma, and W. T. Strunz, Journal of Physics A: Mathematical and Theoretical 54, 035303 (2020).
  * Grossmann (1999) F. Grossmann, Comm. At. Mol. Phys. 34, 141 (1999).
  * Miller (2001) W. H. Miller, J. Phys. Chem. A 105, 2942 (2001).
  * Thoss and Wang (2004) M. Thoss and H. Wang, Ann. Rev. Phys. Chem. 55, 299 (2004).
  * Kay (2005) K. G. Kay, Ann. Rev. Phys. Chem. 56, 255 (2005).
  * Richings  _et al._ (2015) G. W. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt, and B. Lasorne, Int. Rev. in Phys. Chem. 34, 269 (2015).
  * Werther  _et al._ (2021) M. Werther, S. Loho Choudhury, and F. Grossmann, Int. Rev. in Phys. Chem. 40, 81 (2021).
  * Wang  _et al._ (2017) L. Wang, Y. Fujihashi, L. Chen, and Y. Zhao, J. Chem. Phys. 146, 124127 (2017).
  * Werther and Grossmann (2020a) M. Werther and F. Grossmann, Phys. Rev. B 101, 174315 (2020a).
  * Werther and Grossmann (2020b) M. Werther and F. Grossmann, Phys. Rev. A 102, 063710 (2020b).
  * Green and Shalashilin (2019) J. A. Green and D. V. Shalashilin, Phys. Rev. A 100, 013607 (2019).
  * Lee and Fischer (2014) K.-S. Lee and U. R. Fischer, Int. J. Mod. Phys. B 28, 1550021 (2014).
  * Alon  _et al._ (2008) O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Phys. Rev. A 77, 033613 (2008).
  * Simon and Strunz (2014) L. Simon and W. T. Strunz, Phys. Rev. A 89, 052112 (2014).
  * Ray  _et al._ (2016) S. Ray, P. Ostmann, L. Simon, F. Grossmann, and W. T. Strunz, Journal of Physics A: Mathematical and Theoretical 49, 165303 (2016).
  * Viscondi and de Aguiar (2011) T. F. Viscondi and M. A. M. de Aguiar, Journal of Mathematical Physics 52, 052104 (2011).
  * Khripkov  _et al._ (2013) C. Khripkov, D. Cohen, and A. Vardi, Phys. Rev. E 87, 012910 (2013).
  * Rau  _et al._ (2010a) S. Rau, J. Main, and G. Wunner, Phys. Rev. A 82, 023610 (2010a).
  * Rau  _et al._ (2010b) S. Rau, J. Main, H. Cartarius, P. Köberle, and G. Wunner, Phys. Rev. A 82, 023611 (2010b).
  * Polkovnikov  _et al._ (2002) A. Polkovnikov, S. Sachdev, and S. M. Girvin, Phys. Rev. A 66, 053607 (2002).
  * Guaita  _et al._ (2019) T. Guaita, L. Hackl, T. Shi, C. Hubig, E. Demler, and J. I. Cirac, Phys. Rev. B 100, 094529 (2019).
  * Wimberger  _et al._ (2021) S. Wimberger, G. Manganelli, A. Brollo, and L. Salasnich, Phys. Rev. A 103, 023326 (2021).
  * Louisell (1990) W. H. Louisell,  _Quantum Statistical Properties of Radiation_ (Wiley, New York, 1990).
  * Jaksch and Zoller (2005) D. Jaksch and P. Zoller, Annals of Physics 315, 52 (2005), special Issue.
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Greiner  _et al._ (2002) M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature 415 (2002).
  * Kolovsky (2016) A. R. Kolovsky, International Journal of Modern Physics B 30, 1630009 (2016).
  * Perelomov (1971) A. M. Perelomov, Theoretical and Mathematical Physics 6, 156 (1971).
  * Bargmann  _et al._ (1971) V. Bargmann, P. Butera, L. Girardello, and J. R. Klauder, Rep. Math. Phys. 2, 221 (1971).
  * von Neumann (1955) J. von Neumann,  _Mathematical Foundations of Quantum Mechanics_ (Princeton University Press, Princeton, 1955).
  * Boon and Zak (1978) M. Boon and J. Zak, Phys. Rev. B 18, 6744 (1978).
  * Jaksch  _et al._ (1998) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
  * Kramer and Saraceno (1981) P. Kramer and M. Saraceno,  _Geometry of the time-dependent variational principle in quantum mechanics_ (Springer Verlag, Berlin, 1981).
  * Broekhove  _et al._ (1988) J. Broekhove, L. Lathouwers, E. Kesteloot, and P. V. Leuven, Chem. Phys. Lett. 149, 547 (1988).
  * Werther (2020) M. Werther,  _The multi Davydov-Ansatz: Apoptosis of moving Gaussian basis functions with applications to open quantum system dynamics_ , Ph.D. thesis, Technische Universität Dresden (2020).
  * Leggett (2001) A. J. Leggett, Rev. Mod. Phys. 73, 307 (2001).
  * MATLAB (2020) MATLAB,  _9.9.0.1524771 (R2020b) Update 2_ (The MathWorks Inc., Natick, Massachusetts, 2020).
  * Manthe  _et al._ (1992) U. Manthe, H. Meyer, and L. S. Cederbaum, J. Chem. Phys. 97, 3199 (1992).
  * Trotzky  _et al._ (2012) S. Trotzky, Y. Chen, A. Flesch, I. P. McCulloch, U. Schollwöck, J. Eisert, and I. Bloch, Nat. Phys. 8, 325 (2012).
  * Eisert  _et al._ (2010) J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. 82, 277 (2010).
  * Schollwöck (2011) U. Schollwöck, Annals of Physics 326, 96 (2011).