Gaussian time dependent variational principle for the Bose-Hubbard model

TL;DR

This paper extends Gaussian time-dependent variational principle methods to the Bose-Hubbard model, providing improved ground state approximations, excitation spectra, and response functions beyond traditional mean-field theories.

Contribution

It introduces a Gaussian TDVP approach for the Bose-Hubbard model, surpassing mean-field approximations and benchmarking against established methods like DMRG.

Findings

Accurate ground state approximations in 1d, 2d, 3d

Identification of Goldstone, doublon, and continuum excitations

Computed linear response functions relevant to experiments

Abstract

We systematically extend Bogoliubov theory beyond the mean field approximation of the Bose-Hubbard model in the superfluid phase. Our approach is based on the time dependent variational principle applied to the family of all Gaussian states (i.e. Gaussian TDVP). First, we find the best ground state approximation within our variational class using imaginary time evolution in 1d, 2d and 3d. We benchmark our results by comparing to Bogoliubov theory and DMRG in 1d. Second, we compute the approximate 1- and 2-particle excitation spectrum as eigenvalues of the linearized projected equations of motion (linearized TDVP). We find the gapless Goldstone mode, a continuum of 2-particle excitations and a doublon mode. We discuss the relation of the gap between Goldstone mode and 2-particle continuum to the excitation energy of the Higgs mode. Third, we compute linear response functions for perturbations describing density variation and lattice modulation and discuss their relations to experiment. Our methods can be applied to any perturbations that are linear or quadratic in creation/annihilation operators. Finally, we provide a comprehensive overview how our results are related to well-known methods, such as traditional Bogoliubov theory and random phase approximation.

Tables

Table 1. Table 1: Comparison of tangent plane methods. We relate known methods to compute excitation spectra based on choosing a tangent plane of (i) coherent states around | β 0 c ⟩ ket superscript subscript 𝛽 0 c \ket{\beta_{0}^{\mathrm{c}}} , (ii) coherent states around | ψ g ⟩ ket subscript 𝜓 g \ket{\psi_{\mathrm{g}}} and (iii) general Gaussian states around | ψ g ⟩ ket subscript 𝜓 g \ket{\psi_{\mathrm{g}}} and list for which combinations, we find the expected gapless Goldstone mode.

	Linearized TDVP			Projected Hamiltonian
	spectrum of $K^a{}_b$		equivalent method:	spectrum of ${\mathbb{P}}_{|\psi ⟩}\hat{H}{\mathbb{P}}_{|\psi ⟩}$
(i) Coherent states	gapless Goldstone mode ${ℰ}_k^c$	$\Rightarrow$	Bogoliubov theory	gapped 1-particle spectrum
around $|{\beta }_0^c⟩$			(see Appendix A-B)
(ii) Coherent states	gapped 1-particle spectrum ${ℰ}_k^g$	$\Rightarrow$	Iterated Bogoliubov theory	gapped 1-particle spectrum
around $|{\psi }_g⟩$			(see Appendix C)
(iii) Gaussian states	gapless Goldstone mode,	$\Rightarrow$	Random phase approximation	gapped 1-particle spectrum,
around $|{\psi }_g⟩$	gapped 2-particle spectrum		(see Appendix G)	gapped 2-particle spectrum

Equations

$$\hat{H}$$

$$\displaystyle\mathcal{M}=\Big{\{}\ket{\beta,\lambda}=\mathcal{U}(\beta,\lambda)\ket{0}\Big{\}}\,,$$

$$\mathcal{D}(\beta )$$

$$\mathcal{S}(\lambda )$$

$$\displaystyle x^{a}=\big{(}\mathrm{Re}(\beta_{k}),\mathrm{Re}(\lambda_{k,q}),\mathrm{Im}(\beta_{k}),\mathrm{Im}(\lambda_{k,q})\big{)}\,.$$

$${\mathcal{T}}_{|\psi (x)⟩}\mathcal{M}=\text{mbox}span{\left\{\mathcal{U}(x){\hat{b}}_k^{†}|0⟩,\mathcal{U}(x){\hat{b}}_{k-q}^{†}{\hat{b}}_q^{†}|0⟩\right\}}_{k,q}.$$

$$\frac{d}{d\tau }|\psi (\tau )⟩=-{\mathbb{P}}_{|\psi (\tau )⟩}\hat{H}|\psi (\tau )⟩,$$

$$\frac{d}{dt}|\psi (t)⟩={\mathbb{P}}_{|\psi (t)⟩}(-\mathrm{i}\hat{H})|\psi (t)⟩,$$

$${\partial }_V{\mathbb{P}}_{|\psi ⟩}(-\mathrm{i}\hat{H})|\psi ⟩=\left({\partial }_V{\mathbb{P}}_{|\psi ⟩}\right)(-\mathrm{i}\hat{H})|\psi ⟩+{\mathbb{P}}_{|{\psi }_{\mathrm{g}}⟩}(-\mathrm{i}\hat{H})|V⟩,$$

$$K^a{}_b=\frac{\partial }{\partial x^b}\mathrm{Re}⟨W^a|{\mathcal{U}}^{†}(x)\hat{H}\mathcal{U}(x)|0⟩,$$

$$\left\{|W^a⟩\right\}=\left\{{\hat{b}}_k^{†}|0⟩,{\hat{b}}_{k-q}^{†}{\hat{b}}_q^{†}|0⟩,\mathrm{i}{\hat{b}}_k^{†}|0⟩,\mathrm{i}{\hat{b}}_{k-q}^{†}{\hat{b}}_q^{†}|0⟩\right\}.$$

$$K^a{}_b=-\sum _c{\Omega }^{ac}\frac{\partial }{\partial x^c}\frac{\partial }{\partial x^b}E(x),$$

$$\displaystyle\ket{E_{k}}=\mathcal{U}(x^{\mathrm{g}})\,\Big{[}C\,\hat{b}_{k}^{\dagger}+\sum_{q}C_{q}\hat{b}^{\dagger}_{k+q}\hat{b}^{\dagger}_{-q}\Big{]}\ket{0}\,,$$

$$[\hat{H}{]}_{|{\beta }_0^{\mathrm{c}}⟩}=E_{|{\beta }_0^{\mathrm{c}}⟩}-{\Delta }^{\mathrm{c}}+\sum _k{\mathcal{E}}_k^{\mathrm{c}}(\delta {\hat{B}}_k^{\mathrm{c}}{)}^{†}\delta {\hat{B}}_k^{\mathrm{c}},$$

$$\underset{U\to 0}{\lim }\frac{E_{\mathrm{Higgs}}}{U}=\alpha (N,n)$$

$${\hat{H}}_{\lambda }(t)=\hat{H}+\lambda \phi (t)\hat{V},$$

$$\delta V(t)={\frac{d}{d\lambda }⟨{\psi }_{\lambda }(t)⟩\hat{V}{\psi }_{\lambda }(t)}_{\lambda =0}$$

$$Z_V(\omega )=\frac{1}{2}\text{mbox}sign(\omega ){\left|e^a(\omega )d{V}_a\right|}^2\delta (|\omega |),$$

$${\hat{V}}_{1-\mathrm{particle}}^{(k)}$$

$${\hat{V}}_{\mathrm{density}}^{(k)}$$

$${\hat{V}}_{\mathrm{lattice}}^{(k)}$$

$$\hat{H}=\sum _k{\epsilon }_k{\hat{b}}_k^{†}{\hat{b}}_k+\frac{U}{2N}\sum _{k,p,q}{\hat{b}}_{k+q}^{†}{\hat{b}}_{p-q}^{†}{\hat{b}}_k{\hat{b}}_p$$

$${\epsilon }_k=-2\sum _{d=1}^{\dim }\text{mbox}cos\frac{2\pi k_d}{N_d}-\mu ,$$

$$|{\beta }_0⟩=\mathcal{D}({\beta }_0)|0⟩\text{with}\mathcal{D}({\beta }_0)=e^{{\beta }_0{\hat{b}}_0^{†}-{\beta }_0^{\ast }{\hat{b}}_0},$$

$${\beta }_0^{\mathrm{c}}:=\sqrt{-{\epsilon }_{0}N/U}$$

$$E_{|{\beta }_0^{\mathrm{c}}⟩}=⟨{\beta }_0^{\mathrm{c}}|\hat{H}|{\beta }_0^{\mathrm{c}}⟩={\epsilon }_0|{\beta }_0^{\mathrm{c}}{|}^2+\frac{U}{2N}|{\beta }_0^{\mathrm{c}}{|}^4=-\frac{{ϵ}_0^{2}N}{2U},$$

$$[\hat{H}{]}_{|{\beta }_0^{\mathrm{c}}⟩}=E_{|{\beta }_0^{\mathrm{c}}⟩}+\frac{1}{2}\sum _k\left(U_k^{\mathrm{c}}(\delta {\hat{b}}_k^{\mathrm{c}}{)}^{†}\delta {\hat{b}}_k^{\mathrm{c}}+V_k^{\mathrm{c}}\delta {\hat{b}}_k^{\mathrm{c}}\delta {\hat{b}}_{-k}^{\mathrm{c}}+\text{H.c.}\right),$$

$$\mathcal{S}(\lambda )=\exp \left(\frac{1}{2}{\sum }_k({\lambda }_k{\hat{b}}_k^{†}{\hat{b}}_{-k}^{†}-{\lambda }_k^{\ast }{\hat{b}}_k{\hat{b}}_{-k})\right).$$

$$\delta {\hat{B}}_k$$

$$=u_k({\hat{b}}_k-{\beta }_k)-v_k({\hat{b}}_{-k}^{†}-{\beta }_k)$$

$$\tanh 2{\lambda }_k^{\mathrm{c}}=-\frac{V_k^{\mathrm{c}}}{U_k^{\mathrm{c}}},$$

$$[\hat{H}{]}_{|{\beta }_0^{\mathrm{c}}⟩}=E_{|{\beta }_0^{\mathrm{c}}⟩}-{\Delta }^{\mathrm{c}}+\sum _k{\mathcal{E}}_k^{\mathrm{c}}(\delta {\hat{B}}_k^{\mathrm{c}}{)}^{†}\delta {\hat{B}}_k^{\mathrm{c}},$$

Full text

I Introduction

The Bose-Hubbard model provides a theoretical description of interacting cold atoms in optical lattices Bloch et al. (2008), which in the last years have proven to be a promising experimental platform. Its Hamiltonian is given by

[TABLE]

where $\hat{b}^{\dagger}_{i}$ and $\hat{b}_{i}$ are the bosonic creation and annihilation operators for a particle on site $i$ of a square lattice. The model has been analyzed theoretically with a several different methods, ranging from the historical Bogoliubov theory Bogolyubov (1947) to later approaches based on the Gutzwiller ansatz Pekker et al. (2012).

For different choices of the model parameters $U$ and $\mu$, the system exhibits two different phases in the thermodynamic limit: a superfluid phase (small $U$) and a Mott insulator phase (large $U$). One characterization of the superfluid phase is that the $\mathrm{U}(1)$ symmetry generated by the particle number operator $\hat{N}=\sum_{i}\hat{b}_{i}^{\dagger}\hat{b}_{i}$ Sachdev (1999) is spontaneously broken for $N\to\infty$. This leads to both a gapless Goldstone mode and a massive Higgs amplitude mode in the excitation spectrum around the transition. The properties of these have both been described theoretically, e.g., with methods based on the Gutzwiller ansatz Huber et al. (2008, 2007); Bissbort et al. (2014), strong coupling Sengupta and Dupuis (2005); Fitzpatrick and Kennett (2018), the variational cluster approach Knap et al. (2011), the random phase approximation Menotti and Trivedi (2008) or the ladder diagram approximation for the continuum theory Bijlsma and Stoof (1997), and observed in experimental realizations of the model Endres et al. (2012); Bissbort et al. (2011); Schori et al. (2004).

The aim of our paper is to introduce a systematic generalization of the Bogoliubov mean field theory for the superfluid phase. Our method is best described as Gaussian time dependent variational principle (Gaussian TDVP), i.e., we compute system properties from the family of bosonic Gaussian states Weedbrook et al. (2012) given by displaced and squeezed vacua. This is in contrast to Bogoliubov theory which is based on the smaller variational family of coherent states, i.e., just displaced vacua.

Bogoliubov theory describes the model by suitably truncating the Hamiltonian to a quadratic non-interacting mean field Hamiltonian. The minimal energy of this Hamiltonian approximates remarkably well the exact ground state energy. Furthermore, the mean field Hamiltonian can be diagonalized using Bogoliubov transformations and its spectrum describes the dispersion relation of the gapless Goldstone mode of the model. This last step is equivalent to applying the linearized time dependent variational principle to coherent states (coherent TDVP).

Bogoliubov theory, however, also presents several drawbacks. First, the Bogoliubov ground state energy approximation is not variational, i.e., the mean field ground state does not minimize the expectation value with respect to the full Hamiltonian. Second, it does not capture other excitations beyond the Goldstone one, such as the Higgs amplitude mode or bound doublon states. Third, the Goldstone quasiparticles are non-interacting and thus, the decay of quasiparticles excitations can only be studied by re-including the initially discarded Hamiltonian terms as a perturbation Beliaev (1958).

By instead applying linearized TDVP to an extended variational manifold, i.e., the larger class of Gaussian states in place of just coherent states, we overcome all of these drawbacks. First, we compute a variational ground state approximation given by the Gaussian state $\ket{\psi_{\mathrm{g}}}$ with minimal energy expectation value. For this, we use imaginary time evolution and show that $\ket{\psi_{\mathrm{g}}}$ can be efficiently computed in any dimension from two self-consistent equations. Second, our approximate excitation spectrum captures both 1- and 2-particle states, which include the gapless Goldstone mode, a doublon mode and a gapped mode which may be interpreted as a Higgs amplitude mode. The approximate excitation spectrum arises as the eigenvalues of the linearized TDVP equations of motion on the tangent space of Gaussian states. Third, the Gaussian tangent plane naturally captures the interaction of quasi-particle excitations in the 1- and 2-particle sector. This allows us to extract spectral response functions associated to linear and quadratic perturbations and to compute decay and time evolution of excitations.

This paper is structured as follows: In Section II, we introduce our variational manifold and compute the best approximation of the system’s ground state in the superfluid phase, i.e., the Gaussian state with the minimal energy expectation value. In Section III, we study the linearization of the projected real time evolution on such manifold to obtain an expression for the system’s excitation spectrum. In Section IV, we develop a geometric linear response theory consistent with our approximation scheme to capture how linear perturbations couple to different parts of the spectrum. In Section V, we expand on the relationship between our methods and others also based on a Gaussian or coherent state ansatz manifold (specifically Bogoliubov theory) and discuss differences and advantages. Finally, we conclude in Section VI with a comprehensive discussion of our results. In Appendices A-C, we review Bogoliubov theory, its partial equivalence to coherent TDVP and how to make it self-consistent by iteration. In Appendices D-F, we provide further details on the Gaussian ground state approximation, the linearized equations of motion and linear response theory. Finally, in Appendix G we illustrate the equivalence between our Gaussian method and the random phase approximation scheme based on ladder Feynman diagrams.

II Gaussian ground state approximation

As first step of applying our variational methods, we compute the best Gaussian state $\ket{\psi_{\mathrm{g}}}$, i.e., the normalized Gaussian states whose energy expectation value $E_{\ket{\psi_{\mathrm{g}}}}=\bra{\psi_{\mathrm{g}}}\hat{H}\ket{\psi_{\mathrm{g}}}$ on the full Hamiltonian is minimal.

II.1 Variational manifold

We generalize the Bogoliubov theory of the Bose-Hubbard model by extending the variational manifold for the system state to the full manifold $\mathcal{M}$ of bosonic Gaussian states. This is in contrast to regular Bogoliubov theory, where the variation is only done with respect to coherent states. The manifold of Gaussian states can be conveniently parametrized by first squeezing and then displacing the reference vacuum $\ket{0}$, i.e., we consider the variational manifold

[TABLE]

with unitaries $\mathcal{U}(\beta,\lambda)=\mathcal{D}(\beta)\mathcal{S}(\lambda)$ defined by

[TABLE]

where $\hat{b}_{k}=\frac{1}{\sqrt{N}}\sum_{i}e^{-\mathrm{i}kx_{i}}\hat{b}_{i}$ are the momentum space annihilation operators. Here, $\beta_{k}$ is a complex vector and $\lambda_{k,q}$ is a complex matrix invariant under $q\to k-q$. The indices $k$ and $q$ run in the reciprocal lattice. The only redundancy contained in this parametrization is the symmetry of $\lambda$. For a system with $N$ bosonic degrees of freedom, we count $N(N+3)/2$ complex coordinates $(\beta_{k},\lambda_{k,q})$ or $N(N+3)$ real coordinates

[TABLE]

We will use the shorthand notation $\mathcal{U}(x^{\mathrm{g}})$ for the choices of $x^{a}$, such that $\mathcal{U}(x^{\mathrm{g}})\ket{0}=\ket{\psi_{\mathrm{g}}}$.

The manifold is closed under the action of any subgroup generated by any operators that are linear and quadratic in creation/annihilation operators. In particular, this applies to the $\mathrm{U}(1)$ symmetry group generated by the total number operator $\hat{N}=\sum_{i}\hat{b}^{\dagger}\hat{b}_{i}$. Here, for any Gaussian state $\ket{\psi}$ other than the vacuum $\ket{0}$ we find a whole ring of inequivalent states $e^{\mathrm{i}\theta\hat{N}}\ket{\psi}$ with the same energy expectation value. Therefore, we expect this variational manifold to be well suited to capture spontaneously broken $\mathrm{U}(1)$ symmetry phase of the system, i.e., the superfluid phase, and its features, such as the massless Goldstone mode.

While the symmetry of the Bose-Hubbard model is known to be only spontaneously broken in the thermodynamic limit ($N\to\infty$), our ansatz already gives rise to a family of non-symmetric approximate ground states at finite $N$. We furthermore point out that, while finite temperature spontaneous breaking of a continuous symmetry at zero temperature is ruled out in 1D, there might still be quasi-long range order, therefore a broken symmetry ansatz can turn out to be a reasonable choice also in 1D.

The manifold contains a submanifold of states which are translationally invariant, namely the set of states $\ket{\psi(\beta_{k},\lambda_{k,q})}$ with $\beta_{k}=\delta_{k,0}\beta_{0}$ and $\lambda_{k,q}=\delta_{k,0}\lambda_{0,q}$. For the ground state search, it is sufficient to restrict ourselves to this submanifold as we expect the ground state to preserve the translational symmetry of the problem. For the study of excitations around the translationally invariant ground state, we will then use the full manifold in order to capture also excitations with non-zero momentum.

The tangent space $\mathcal{T}_{\ket{\psi}}$ of the variational manifold at the point $\ket{\psi}$ is naturally spanned by the states with 1- and 2-particle excitations, i.e.,

[TABLE]

Put differently, the variational class of all Gaussian states captures accurately the 1- and 2-particle quasiparticle excitation sector of our model.

II.2 Imaginary time evolution

The first step of our procedure to exploit the given choice of variational manifold is to find within it the best approximation of the ground state, that is the state with the lowest energy expectation value. One strategy to do this is to consider the projected imaginary time evolution. This is the solution of the evolution equation

[TABLE]

where $\mathbb{P}_{\ket{\psi(\tau)}}$ is the orthogonal projector onto the tangent space to the manifold at $\ket{\psi(\tau)}$. This projection ensures that the solution will be contained in the ansatz manifold at all times $\tau$. This evolution converges from a random initial state to a local minimum of the energy expectation value function and can be shown to be equivalent to a gradient descent method. For Gaussian states, we find simple equations for the stationary point of this evolution, i.e., the state $\ket{\psi_{\mathrm{g}}}\in\mathcal{M}$ such that $-\mathbb{P}_{\ket{\psi_{\mathrm{g}}}}\hat{H}\ket{\psi_{\mathrm{g}}}=0$, and see that they only admit a single solution up to the redundancy generated by $e^{\mathrm{i}\alpha\hat{N}}$.

This solution $\ket{\psi_{\mathrm{g}}}$ can be characterized analytically, independently of the system size or dimensionality, in terms of two parameters $A$ and $B$, which can be efficiently computed numerically as the fixed point of two coupled self-consistent equations. For more details on this calculation and on how to parametrize the approximate ground state see Appendix D.

II.3 Ground state properties

Having obtained an analytical expression for the approximate ground state, it is then possible to calculate the predictions of our model for ground state properties such as the energy and particle densities. The quality of our method can be benchmarked by comparing these quantities with the ones obtained through other methods, such as Bogoliubov theory or, at least in one dimension, with a numerical DMRG White (1992) calculation (see Figure 1). Our variational energy $E_{\ket{\psi_{\mathrm{g}}}}=\bra{\psi_{\mathrm{g}}}\hat{H}\ket{\psi_{\mathrm{g}}}$ is higher than the DMRG one, as expected, but lower than the one obtained by other variational choices, such as the coherent state $\ket{\beta^{\mathrm{c}}_{0}}$ with minimal energy $E_{\ket{\beta^{\mathrm{c}}_{0}}}$. The energy obtained as the ground state energy of the Bogoliubov mean field Hamiltonian is generally lower than ours and remarkably close to the DMRG result. However, it is important to emphasize that this energy $E_{\mathrm{Bogoliubov}}$ is not variational as it is computed with respect to the truncated mean field Hamiltonian, which actually does not admit a well defined ground state in the zero momentum mode. More precisely the state minimizing the mean field energy is infinitely squeezed, which would lead to a diverging energy expectation value with respect to the full Bose-Hubbard Hamiltonian.

Thus, the Gaussian variational family provides a consistent class to approximate the ground state of the Bose-Hubbard model in the superfluid phase, even though its ground state energy estimate is worse than the one obtained from Bogoliubov theory. However, the strength of our extended variational family lies in its prediction of quasiparticle excitations and their properties, such as life time and linear response.

III Quasi particle excitations

We can derive an approximate excitation spectrum from the perspective of our Gaussian variational manifold by looking at real time evolution of the system projected onto the manifold. As our variational class generalizes the set of coherent state manifold used in standard Bogoliubov theory, we will be able to capture higher excitation modes of the model.

The projected real time evolution is computed as prescribed by the time dependent variational principle (TDVP) Haegeman et al. (2011); Dirac (1930). Such principle can be formulated as stating that the real time evolution projected on the manifold of Gaussian states takes the form

[TABLE]

and generates a Hamiltonian time evolution flow $\Phi_{t}:\mathcal{M}\rightarrow\mathcal{M}$ that, in a neighbourhood of the stationary state $\ket{\psi_{\mathrm{g}}}$, reduces to a sum of phase rotations. From the perspective of our variational manifold, the frequencies of these rotations provide a natural approximation of the lowest excitation energies.

III.1 Linearized TDVP

We calculate the excitation energies as shown in Figure 2, as the eigenvalues of the linearized equations of motion, that can be understood as a generalization of the well-known Gross-Pitaevskii equation Pitaevskii (1961); Gross (1961).

In particular, we consider the linearization of this projected real time evolution around the stationary point $\ket{\psi_{\mathrm{g}}}$, that defines a linear map $K$ at the tangent space of $\ket{\psi_{\mathrm{g}}}$, mapping a tangent vector $\ket{V}$ onto

[TABLE]

where $\partial_{V}$ indicates the directional derivative on $\mathcal{M}$ evaluated at $\ket{\psi_{\mathrm{g}}}$ in the direction indicated by $\ket{V}$.

Using the fact that $\ket{\psi_{\mathrm{g}}}$ is a stationary point of the evolution, i.e., $\mathbb{P}_{\ket{\psi_{\mathrm{g}}}}(-\mathrm{i}\hat{H})\ket{\psi_{\mathrm{g}}}=0$, we can write the map $K$ as a real $N(N+3)\times N(N+3)$ matrix, referring to the tangent plane basis $\{\mathcal{U}(x^{\mathrm{g}})\ket{W^{a}}\}$ (see Appendix E),

[TABLE]

where $x^{a}$ refers to the real coordinates (5) of the manifold and the orthonormal frame $\ket{W^{a}}$ is given by

[TABLE]

To give this matrix representation of $K$ we are here employing a real formalism in which we use a real parametrization and we consider the tangent plane as real vector space, i.e., we consider the vectors $\ket{W}$ and $\ket{W^{\prime}}=\mathrm{i}\!\ket{W}$ to be linearly independent (and orthogonal with respect to the real inner product $\mathrm{Re}\braket{W}{W^{\prime}}$). This is naturally dictated by the fact that the linearized equations of motion from (9) do not commute with complex multiplication, i.e., $K^{a}{}_{b}$ is not complex-linear, as elaborated in Appendix E. Another formal expression for the matrix $K$ is

[TABLE]

where $E(x)$ is the energy expectation value of the state $\ket{\psi(x)}$, and the matrix $\Omega$ is the antisymmetric symplectic form defined in (108).

The evaluation of this matrix (10) reduces to calculating expectation values using Wick’s theorem and taking derivatives, therefore it can be calculated analytically in terms of the ground state parameters obtained in the previous section. More details on the form of $K$ can be found in Appendix E. $K$ is a symplectic matrix whose eigenvalues come in complex conjugate pairs $\pm\mathrm{i}\omega$. The values of $\omega$ are our estimates of the excitation energies of the model.

Due to the translational invariance of $\hat{H}$ the matrix $K$ is block diagonal, with each block acting on the span of tangent vectors with fixed total momentum, which we labeled by $k$ in in (92). The approximate excitation energies $\omega$ can therefore also be labeled by the total momentum $k$ their respective eigenvector. The size of each block grows linearly in $N$, and therefore in the thermodynamic limit $N\to\infty$ there is an infinity of eigenvalues $\omega_{k}$ for each $k$, which can arrange themselves in a continuum plus possibly some discrete excitations that represent bound states.

III.2 Excitation spectrum

In Figure 2, we show the dispersion relations obtained by diagonalizing the matrix $K$ numerically. For momenta close to zero, we always find a gapless isolated mode that agrees well with the Bogoliubov dispersion relation $\mathcal{E}_{k}$. However, we also find a continuum of states that have energies above this Goldstone mode and that always shows a gap around $k=0$. Finally, for certain parameter choices, e.g., for strong interactions, our spectrum also contains another isolated state above the continuum, which can be interpreted as a doublon state.

We point out that the fact that our method gives a gapless mode was to be expected. Indeed, $\hat{N}$ commutes with the Hamiltonian and the vector $\hat{N}\ket{\psi}$ is part of the tangent plane for all $\ket{\psi}\in\mathcal{M}$, because $\hat{N}$ is quadratic in the bosonic creation and annihilation operators. Therefore there exists a direction in the manifold along which the energy is constant. In this direction, the Hessian $\frac{\partial}{\partial x^{a}}\,\frac{\partial}{\partial x^{b}}E(x)$ has a vanishing eigenvalue and thus, because of equation (12), also $K$ does.

Our method captures the tangent plane generated by displacements and squeezing, i.e., it is spanned by 1- and 2-particle excitations. A generic eigenvector $\ket{E_{k}}$ of $K^{a}{}_{b}$ with momentum $k$ is

[TABLE]

where $C,C_{q}\in\mathbb{C}$. We should therefore compare our results with the 1- and 2-particle excitation spectrum obtained from Bogoliubov theory.

Traditional Bogoliubov theory constructs the excitation spectrum from the 1-particle dispersion relation $\mathcal{E}_{k}$ (see (34)) of the mean field Hamiltonian

[TABLE]

as reviewed in A. The dispersion relation $\mathcal{E}_{k}$ is independent of the interaction strength $U$ and becomes exact in the limit $U\to 0^{+}$. General eigenstates of $[\hat{H}]_{\ket{\beta^{\mathrm{c}}_{0}}}$ consist of non-interacting excitations created by $(\delta\hat{B}^{\mathrm{c}}_{k})^{\dagger}$. A general 2-particle excitation with momentum $k$ is therefore given by $(\delta\hat{B}_{k+q}^{\mathrm{c}})^{\dagger}(\delta\hat{B}^{\mathrm{c}}_{k-q})^{\dagger}\ket{\beta^{\mathrm{c}}_{0}}$ and has energy $\mathcal{E}_{k+q}+\mathcal{E}_{k-q}$.

Because of the gapless nature of the 1-particle Bogoliubov dispersion relation, the continuum of non-interacting 2-particle excitations is never separated in energy from the 1-particle dispersion, as seen in Figure 3 (right) and discussed in Appendix A. The gap between the isolated bound state (Goldstone mode) and the continuum of higher excitations is therefore a new feature of Gaussian TDVP due to the fact that it implements the interaction within the 1- and 2-particle sectors.

While the Goldstone mode continues to be well-described by the Bogoliubov dispersion relation $\mathcal{E}_{k}$, the spectrum of 2-particle excitations from Gaussian TDVP starts to divert as we increase $U$. In particular, we see that for sufficiently large $U$ both, the Goldstone mode and the doublon mode are completely separate from the continuum.

In figure 3, we compare Gaussian TDVP and Bogoliubov theory in the regime where the Goldstone mode partially intersects with the continuum. We observe that this intersection appears for small $U$ in the Gaussian TDVP results only in those regimes where also the in Bogoliubov theory the 1-particle mode lies partially above the bottom of the many particle continuum. This phenomenon occurs for choices of $\mu$ and system dimension $\dim$ such that the dispersion relation $\mathcal{E}_{k}$ is not convex, i.e., there exist $q,k$, such that $\mathcal{E}_{k}+\mathcal{E}_{q}<\mathcal{E}_{q+k}$. In Appendix A, we show that this can only happen for $\mu<6-2\dim$. The Gaussian TDVP continuum (light orange) agrees well with the Bogoliubov 2-particle spectrum (light red), where it intersects with the Goldstone mode, i.e., roughly for $2\pi k/N\in(\pi/4,3\pi/4)$. Outside of this region, the two disagree: While the Gaussian TDVP gives rise to a finite gap between continuum and isolated Goldstone mode, 2-particle continuum and 1-particle dispersion relation necessarily touch for the non-interacting mean field Hamiltonian from Bogoliubov theory.

In this regime of intersecting continuum and Goldstone mode, Gaussian TDVP can (at least partially) describe the decay of the Goldstone mode into the continuum of modes with a higher number of excitations. This phenomenon is known as Beliaev damping Beliaev (1958) and is not captured by the standard Bogoliubov theory, but so far has been typically obtained from perturbative expansions by re-including higher order terms of the Hamiltonian. We will investigate this decay of 1-particle excitations into the continuum in section IV.2.

III.3 Higgs mode

Another suggestive observation can be made on the physical interpretation of the gapped mode at the bottom of the continuum. A possible interpretation is that it is a remnant of what, near the superfluid to Mott insulator transition, becomes known as the Higgs mode. It corresponds to oscillations of the amplitude of the order parameter $\braket{b_{0}}$ (while the Goldstone mode is interpreted as oscillations of the order parameter phase) and it has been observed experimentally by coupling to it through modulation of the tunneling amplitude Endres et al. (2012).

The prediction of our model for the continuum gap, which because of this possible interpretation we will label as $E_{\mathrm{Higgs}}$, can be studied numerically through the diagonalization of the matrix $K$ described in the previous paragraphs. At fixed non-zero $U$, $E_{\mathrm{Higgs}}$ converges to a finite non-zero value in the thermodynamic limit.

We are also able to give an analytical asymptotic result for the limit in which $U\to 0$ while $\mu$ varies so as to keep a constant particle number density of the ground state $n=\braket{\hat{N}}/N$ (see Appendix E). In this limit, we have that the gap goes to zero linearly in the interaction strength $U$, namely

[TABLE]

Note that it is instrumental that we took here first the limit $U\to 0$, before studying the large $N$ asymptotics.

In Figure 4 one can see the numerical results for the behaviour of the ratio between the Higgs gap and $U$ and notice how it indeed approaches a constant asymptotic value for small $U$. In the large $U$ region, it has instead an unexpected divergent behaviour (the gap should close at the SF/MI transition Huber et al. (2007)), however this can be understood as a breaking down of our model at the transition where Gaussian states are no longer a good description of the system’s ground state. It is instead interesting to see how the constant small $U$ behaviour matches the experimentally measured value of the Higgs mode gap Endres et al. (2012) even better than the previous theoretical results obtained with Gutzwiller theory.

IV Linear Response

We use our variational manifold and the real time evolution projected onto it to study the response of the system to small perturbations. This is significant as it provides possible connections to actual experiments, where certain system responses can be probed and measured.

IV.1 Spectral Functions

We model an external perturbation by considering the time dependent perturbed Hamiltonian

[TABLE]

where $\hat{H}$ is the unperturbed Bose-Hubbard Hamiltonian (1), $\varphi(t)$ is a classical external field that couples to the system through the Hermitian operator $\hat{V}$ and $\lambda$ is a real parameter. We shall then consider the projected real time evolution $\ket{\psi_{\lambda}(t)}$ of the system under such perturbed Hamiltonian and evaluate its response in terms of the expectation value of the same coupling operator $\hat{V}$. In particular, we consider this response in the limit of small perturbations, i.e., we compute quantities only up to first order in the parameter $\lambda$. Thus, we consider the response

[TABLE]

to the perturbation $\hat{V}$.

As discussed in Appendix F, the Fourier space response $\delta V(\omega)$, calculated on the variational manifold as explained above, takes the form $\delta V(\omega)=\tilde{\varphi}(\omega)\,\chi(\omega)$, where $\tilde{\varphi}$ is the Fourier transform of the perturbing field $\varphi(t)$ and $Z_{V}(\omega)\equiv-\pi\mbox{Im}\chi(\omega)$ is the response function of the system with respect to the perturbation $\hat{V}$. Such response functions are expressed in terms of the spectral decomposition of the linearized real time evolution $K$ defined in Section III as

[TABLE]

where $dV$ is the gradient differential form of the real valued function on the manifold $\braket{\psi(\beta,h)}{\hat{V}}{\psi(\beta,h)}$, $e(\omega)$ are the eigenvectors of $K$ (defined in Appendix E) with eigenvalue $\omega$ and $\delta(|\omega|)$ is a normalization of the eigenvectors defined in equation (121).

For the Bose-Hubbard model, we consider the following types of perturbations:

[TABLE]

In (19), we use a linear operator to create a single particle perturbation of momentum $k$. The other two perturbations are quadratic in creation and annihilation operators, such that the excitation consists in general of both single- and 2-particle excitations. In (20), we consider a spatial density variation by modulating the chemical potential with momentum $k$, which couples directly to the local particle density. In (21), we consider a spatial modulation of momentum $k$ of the hopping constant. This can be achieved through a modulation of the lattice depth Kollath et al. (2006). Such perturbation naturally couples to the kinetic energy operator.

The different response functions $Z_{k}(\omega)$, obtained by evaluating (18) for different types of perturbation operators of momentum $k$ and at energy $\omega$, give us an indication of how strongly each type of perturbation couples to different regions of the spectrum.

A first observation we can make is on the behaviour of the isolated Goldstone mode in those situations when it merges with the continuum part of the spectrum. In Figure 7, we see how the isolated peak of the response function broadens into a wider feature inside the continuum. This indicates how, even when the Goldstone mode is not an isolated eigenstate, it still survives as a finite lifetime excitation of the system.

We can then also compare how the different perturbations considered couple to the system. In Figure 6, we see how the perturbation that couples the strongest to the continuum modes at $k=0$ is the lattice modulation operator. Although it has to be mentioned that the definition of the normalizations of the perturbations (20-21) is not free of some arbitrariness, the large difference in these coupling strengths provides a further element of support to the identification of the lower continuum modes as the Higgs excitation. Indeed it is known that the Higgs mode should be excited most easily through perturbations of the kinetic energy term of the Hamiltonian, while the Goldstone mode through perturbations in the particle density Huber et al. (2007).

IV.2 Real time evolution

The analysis of the response function can also give indications on the real time evolution of perturbations of the system. Indeed, we can interpret the operator $\hat{V}$ as creating a perturbation described by the tangent vector $\ket{\delta\psi_{\mathrm{g}}(0)}=\mathbb{P}_{\ket{\psi_{\mathrm{g}}}}(-\mathrm{i}\hat{V})\ket{\psi_{\mathrm{g}}}$ at $t=0$, which is equivalent to giving the system a kick by choosing $\varphi(t)=\delta(t)$. The evolution in the tangent plane of this perturbation vector is then given by $d\Phi_{t}$ at $\ket{\psi_{\mathrm{g}}}$, i.e., the push-forward of the real time evolution flow $\Phi_{t}$ around the stationary point. $d\Phi_{t}$ is a linear map on the tangent space at $\ket{\psi_{\mathrm{g}}}$, explicitly it is given by the matrix $e^{Kt}$ with respect to the basis (6).

We consider an initial perturbation $\ket{\delta\psi_{\mathrm{g}}(0)}$ in the 1-particle sector of the tangent plane, i.e., we require $\ket{\delta\psi_{\mathrm{g}}(0)}$ to be only spanned by 1-particle states in basis (6). In particular, this is accomplished by the perturbation created by $\hat{V}^{(k)}_{1-\mathrm{particle}}$. The time evolution of $\ket{\delta\psi_{\mathrm{g}}(0)}$ under the map $d\Phi_{t}$ can then show different characteristic behaviour. If the perturbation is created at a momentum value where there exists an isolated Goldstone state with a strong coupling to the 1-particle sector (quantified by the 1-particle response function), the perturbation will persist indefinitely. If, instead, the excitation has a momentum at which the 1-particle perturbation couples sufficiently strongly to the continuum, a part of it will decay into the continuum modes, disappearing in a time proportional to the inverse of the width of the response function. Finally, if there is no isolated Goldstone state at the chosen momentum, but only the continuum, the perturbation will have a finite lifetime and decay completely into continuum excitations.

In Figure 8, we show the overlap of the time evolved perturbation with the 1-particle sector of the tangent plane, for different total momenta of the initial perturbation. A perturbation with momentum $k$ corresponding to an isolated Goldstone state will maintain a large overlap with the 1-particle sector. For perturbations with momentum $k$ closer to the region where the Goldstone mode merges with the continuum, a larger part of the overlap with the 1-particle sector will decay in time. Finally, if the perturbation has a total momentum $k$, for which no isolated Goldstone state exists in the spectrum, the single particle overlap will decay completely to zero after a finite lifetime. Such decay behavior is similar to what can also be seen in quantum optical systems coupled to unconventional photon baths González-Tudela and Cirac (2017).

This behaviour of the evolution of perturbations can be interpreted as a remnant in lattice systems of what in continuum Bose-Einstein condensates (BEC) is known as the Beliaev damping of excitations, i.e., the decay of 1-particle excitations into the continuum of many particle excitations due to scattering interactions. Our variational scheme successfully captures at least part of this behaviour, namely the one associated to the 1- and 2-particle sector that are fully included in our tangent space. This is in contrast to the traditional Bogoliubov theory that is restricted to the non-interacting 1-particle sector. In particular, standard Bogoliubov theory cannot describe the interaction with the continuum consistently, which can only be incorporated by re-including the previously neglected terms as perturbations Beliaev (1958).

V Relations between methods

Our study is based on the time dependent variational principle (TDVP), where we project the equations of motion on a given variational class and linearize them around the stationary state that provides the best approximation of the ground state. While we focused on the class of all Gaussian states, the method can be applied to any suitable family of states, so it is natural to compare the results between different variational classes. In the context of Bogoliubov theory, it is natural to compare our larger manifold of all bosonic Gaussian states $\mathcal{D}(\beta)\mathcal{S}(\lambda)\ket{0}$ with the smaller sub manifold consisting only of coherent (or displaced) states $\mathcal{D}(\beta)\ket{0}$. Table 1 summarizes the different methods.

(i) Coherent TDVP around $\ket{\beta_{0}^{\mathrm{c}}}$. If we apply linearized TDVP to the manifold of coherent states, we obtain the same excitation spectrum as the single-particle spectrum of Bogoliubov mean field theory (see Appendix B). The latter is defined by taking the full Hamiltonian and using the commutation relations to normal order the creation and annihilation operators with respect to the coherent state $\ket{\beta_{0}^{\mathrm{c}}}$ that minimizes the energy on the coherent state manifold. At this point, we can truncate at quadratic order to obtain the mean field Hamiltonian $[\hat{H}]_{\ket{\beta_{0}^{\mathrm{c}}}}$ and use a Bogoliubov transformation to compute its excitation spectrum $\mathcal{E}_{k}$. However, we should point out that $[\hat{H}]_{\ket{\beta_{0}^{\mathrm{c}}}}$ contains more information than the linearized TDVP, as it gives us a Hamiltonian operator whose minimal energy $E_{\mathrm{Bogoliubov}}$ is a better estimate of the system’s ground state energy than just $E_{\ket{\beta_{0}^{\mathrm{c}}}}$. On the other hand, this energy is not variational, i.e., it cannot be expressed as the expectation value of an ansatz state on the full system Hamiltonian. Furthermore, the truncation of $[\hat{H}]_{\ket{\beta_{0}^{\mathrm{c}}}}$ is not self-consistent, because $\ket{\beta_{0}^{\mathrm{c}}}$ is not its ground state.

(ii) Coherent TDVP around $\ket{\psi_{0}^{\mathrm{c}}}$. After finding the best Gaussian ground state approximation $\ket{\psi_{\mathrm{g}}}$, we can linearize the equations of motion restricted to the space of displacements. This is equivalent to iterating traditional Bogoliubov theory as reviewed in appendix C, where we find the self-consistent mean field Hamiltonian $[\hat{H}]_{\ket{\psi_{\mathrm{g}}}}$, whose ground state is again $\ket{\psi_{\mathrm{g}}}$. The 1-particle spectrum $\mathcal{E}_{k}^{\mathrm{g}}$ from this Hamiltonian is gapped and consequently not a good approximation to the Goldstone mode. However, we can use $\mathcal{E}_{k}^{\mathrm{g}}$ to construct the 2-particle continuum of the quadratic Hamiltonian $[\hat{H}]_{\ket{\psi_{\mathrm{g}}}}$. Interestingly, the resulting 2-particle spectrum provides a good approximation to the continuum with Gaussian TDVP (see figure 11). In this way, we can understand the gap $\mathcal{E}^{\mathrm{g}}_{0}$ as already encoding the interaction energies between two particle excitations that is required to approximate the interacting 2-particle spectrum.

(iii) Gaussian TDVP around $\ket{\beta_{0}^{\mathrm{c}}}$. In order to obtain a self-consistent ground state, we enlarge the manifold of states and introduce general Gaussian states, which also allow for squeezing. Indeed, the Gaussian state of minimal energy $\ket{\psi_{\mathrm{g}}}$ can also be identified as the state that fulfills the property of being the ground state of the corresponding mean field Hamiltonian $[\hat{H}]_{\ket{\psi_{\mathrm{g}}}}$, i.e., the quadratic truncation of the full Hamiltonian when normal-ordered with respect to $\ket{\psi_{\mathrm{g}}}$. If we apply linearized TDVP to the extended manifold of Gaussian states, we obtain the spectrum object of this paper, which naturally contains both 1- and 2-particle excitations (see Section III). The Gaussian TDVP spectrum can be equivalently obtained using random phase approximation (see Appendix G). More precisely, the Gaussian state $\ket{\psi_{\mathrm{g}}}$ can be used as the reference vacuum when expanding Green’s functions in terms of Feynman diagrams. We can consistently resum all ladder diagrams to obtain an approximate excitation spectrum, that agrees with the one found through Gaussian TDVP.

Projected Hamiltonian. Finally, there is a well-known alternative Haegeman et al. (2013); Vanderstraeten et al. (2019); Shi et al. (2018) to compute excitation spectra from a tangent plane based on the projected Hamiltonian. Instead of linearizing the equations of motion, we can directly take the tangent plane as variational ansatz for eigenstates by projecting the full Hamiltonian onto it, i.e., $H_{\mathbb{P}}=\mathbb{P}_{\ket{\psi}}\hat{H}\mathbb{P}_{\ket{\psi}}$, and then computing its spectrum. The eigenstates $\ket{E_{i}}$ with energy $E_{i}$ of the projected Hamiltonian $H_{\mathbb{P}}$ are manifestly variational, i.e., their expectation value with respect to the full Hamiltonian is equal to $E_{i}$ and there exists a true eigenstate of the full Hamiltonian with smaller energy. This is not necessarily the case for the eigenvectors of $K^{a}{}_{b}$ in the linearized TDVP. In Haegeman et al. (2013), it has been further pointed out that–in contrast to the projected Hamiltonian method–the linearized TDVP may incorrectly predict massless excitation modes. This occurs whenever the approximate ground state within the chosen variational family spontaneously breaks a symmetry which is not spontaneously broken in the exact ground state. In the case of the Bose-Hubbard model, this is actually a desirable feature: while the true ground state only breaks the $\mathbb{U}(1)$ symmetry in the limit $N\to\infty$, the family of Gaussian states already breaks this symmetry for finite $N$ and is thus well-suited to study the superfluid phase in the thermodynamic limit.

VI Discussion and outlook

Our Gaussian TDVP method naturally generalizes Bogoliubov theory to describe the superfluid phase of the Bose-Hubbard model. The presented methods provide systematic framework to compute (a) approximate ground state energies, (b) excitation spectra and (c) linear response functions for general variational families.

(a) Our variational ansatz of all Gaussian states provided a good variational approximation of the degenerate set of ground states. While its ground energy prediction is generally worse than the ground state energy of the traditional Bogoliubov mean field Hamiltonian, our method has the advantage of being manifestly variational, i.e., we find a concrete state whose energy expectation value with respect to the full Hamiltonian provides a rigorous upper bound to the true ground state energy.

(b) The linearized equations of motion projected onto the full Gaussian tangent plane allowed us to capture interactions between 1-particle excitations. We find a gapless mode that can be identified with the Goldstone mode that is well-approximated by traditional Bogoliubov theory. We point out that, compared to other variational approaches such as Knap et al. (2011), our variational family is chosen so as to yield a Goldstone mode that is exactly gapless. We also find a gap between Goldstone mode and 2-particle continuum that is not captured by traditional Bogoliubov theory’s non-interacting $n$-particle spectrum. We argued that the lowest band of the 2-particle continuum can be identified as the Higgs mode. Finally, for some parameter choices, we found an isolated doublon state.

(c) Using eigenvalues and eigenvectors of the linearized equations of motion, we could compute spectral functions for linear perturbations that are generated by arbitrary linear and quadratic operators. Here, we observed how these different perturbations coupled to isolated bound states, e.g., the Goldstone and the doublon mode. Finally, we described the time evolution of the Goldstone mode, in particular its (partial) decay depending on its interplay with the 2-particle continuum.

An interesting feature of the TDVP is the interplay and relations to other methods, such as Bogoliubov theory and random phase approximation. For the computation of excitation spectra, we find that these standard methods give the same results as the TDVP computation for the correct choice of variational manifold. Another interesting observation is the fact that the 2-particle continuum of the Gaussian TDVP is well-approximated by the free 2-particle spectrum from coherent TDVP (or equivalently: iterated Bogoliubov theory) around the best Gaussian state $\ket{\psi_{\mathrm{g}}}$. This suggests that the 1- and 2-particle spectrum can be approximated by a hybrid approach with two different quadratic Hamiltonians, i.e., the Bogoliubov Hamiltonian $[\hat{H}]_{\ket{\beta_{0}^{\mathrm{c}}}}$ and the iterated Bogoliubov Hamiltonian $[\hat{H}]_{\ket{\psi_{\mathrm{g}}}}$. While the former describes the massless Goldstone mode, we can use the latter’s free 2-particle spectrum to approximate the interacting 2-particle spectrum of the full Hamiltonian.

The one of our predictions that calls most for further inquiry is the gapped $2$-particle continuum above the Goldstone mode. It will be interesting to further explore with other methods whether the identification of the lowest continuum mode as the Higgs mode is correct and whether the gap that separates it from the Goldstone mode survives, once one also considers excitations of three or more particles.

The presented scheme is self-consistent and requires no other assumptions than the choice of variational manifold. In particular, it can be easily applied to other variational families, such as non-Gaussian states Shi et al. (2018) or Gutzwiller states. To study the Bose-Hubbard model also in the Mott phase, we expect that a variational manifold that combines both, Gutzwiller states and Gaussian transformations, i.e., the family of states resulting from applying a Gaussian unitary to a Gutzwiller product state, is particularly promising and will be analyzed in future work.

Acknowledgements.

We thank Abhay Ashtekar, Bertrand Halperin, Andy Martin, Marcos Rigol, Richard Schmidt and Yao Wang for inspiring discussions. TG, LH and IC thank Harvard University for the hospitality during several visits. TG, LH, CH and IC are supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2111 – 39081486. LH is funded by the the Max Planck Harvard Research Center for Quantum Optics. TS acknowledges the Thousand-Youth-Talent Program of China. CH and IC acknowledge funding through ERC Grant QUENOCOBA, ERC-2016-ADG (Grant no.742102). ED acknowledges funding through Harvard-MIT CUA, AFOSR-MURI: Photonic Quantum Matter (award FA95501610323) and DARPA DRINQS program (award D18AC00014).

Appendix A Review of Bogoliubov theory

In the main part of this paper, we generalize the well-known Bogoliubov theory of the Bose-Hubbard model by extending the underlying variational family from coherent to squeezed coherent states. To allow for a fair comparison, we now review the three steps involved in traditional Bogoliubov theory.

For later computations, it is useful to write the Bose-Hubbard Hamiltonian (1) in momentum space

[TABLE]

where we defined $\hat{b}_{k}=\frac{1}{\sqrt{N}}\sum_{i}e^{-\mathrm{i}kx_{i}}\hat{b}_{i}$ on the reciprocal lattice and introduced the non-interacting dispersion relation

[TABLE]

where $N_{d}$ refers to the number of lattice sites in the $d$-th direction, such that $N=\prod_{d}N_{d}$.

Step 1 (coherent variation). Bogoliubov theory approximates the ground state within the class of translationally invariant coherent states, i.e., the states

[TABLE]

satisfying $\bra{\beta_{0}}b_{k}\ket{\beta_{0}}=\beta_{0}\delta_{0,k}$. Within this class, the average energy value is minimized for $|\beta_{0}|$ equal to

[TABLE]

leading to the expectation value

[TABLE]

provided that $\varepsilon_{0}<0$. There is a larger set of solutions given by $\beta_{0}=\beta_{0}^{\mathrm{c}}e^{\mathrm{i}\varphi}$ associated to the spontaneously broken $\mathrm{U}(1)$ symmetry generated by $\hat{N}=\sum_{k}\hat{b}_{k}^{\dagger}\hat{b}_{k}$.

Step 2 (mean field Hamiltonian). We can use $\ket{\beta^{\mathrm{c}}_{0}}$ to define the mean field Hamiltonian

[TABLE]

where $\delta\hat{b}_{k}^{\mathrm{c}}=\hat{b}_{k}-\delta_{k,0}\beta^{\mathrm{c}}_{0}$, $U_{k}^{\mathrm{c}}=\varepsilon_{k}+\frac{2U}{N}\,|\beta^{\mathrm{c}}_{0}|^{2}=\varepsilon_{k}-2\varepsilon_{0}$ and $V_{k}^{\mathrm{c}}=\frac{U}{N}\,(\beta^{\mathrm{c}}_{0})^{2}=-\varepsilon_{0}$. Here, we define $[\hat{H}]_{\ket{\psi}}$ to be the quadratic Hamiltonian resulting from the quadratic truncation of $\hat{H}$ written as normal ordered polynomial in creation and annihilation operators $\delta\hat{b}_{k}^{\dagger}$ and $\delta\hat{b}_{k}$ associated to the Gaussian state $\ket{\psi}$, i.e., in our case, $\delta\hat{b}_{k}\ket{\beta^{\mathrm{c}}_{0}}=0$.

Step 3 (squeezed ground state). The mean field Hamiltonian $[\hat{H}]_{\ket{\beta^{\mathrm{c}}_{0}}}$ is quadratic, implying that we can diagonalize it by applying the Bogoliubov transformation

[TABLE]

We perform the transformation by expressing $\hat{b}_{k}$ in terms of new creation and annihilation operators

[TABLE]

The diagonalization is accomplished by $\beta^{\mathrm{c}}_{0}$, $u_{k}^{\mathrm{c}}=\cosh{\lambda_{k}^{\mathrm{c}}}$ and $v_{k}^{\mathrm{c}}=\sinh{\lambda_{k}^{\mathrm{c}}}$, where $\lambda^{\mathrm{c}}_{k}$ is given by

[TABLE]

such that the Hamiltonian takes the form

[TABLE]

where the energy shift $\Delta^{\mathrm{c}}$ and the excitations $\mathcal{E}_{k}^{\mathrm{c}}$ are

[TABLE]

The formal ground state of this Hamiltonian is given by a state $\ket{\beta_{0}^{\mathrm{c}},\lambda^{\mathrm{c}}}=\mathcal{D}(\beta_{0}^{\mathrm{c}})\mathcal{S}(\lambda^{\mathrm{c}})\ket{0}$, such that $\delta\hat{B}^{\mathrm{c}}_{k}\ket{\beta_{0}^{\mathrm{c}},\lambda^{\mathrm{c}}}=0$. The Bogoliubov spectrum is gapless, i.e., we have $\mathcal{E}^{\mathrm{c}}_{0}=0$.

The full spectrum of the mean field Hamiltonian $[\hat{H}]_{\ket{\beta^{\mathrm{c}}_{0}}}$ is constructed from the 1-particle Bogoliubov spectrum $\mathcal{E}^{\mathrm{c}}_{k}$, whose excitations are non-interacting. Depending on the shape of the function $\mathcal{E}^{\mathrm{c}}_{k}$, i.e., if there exist $k,q$ with $\mathcal{E}^{\mathrm{c}}_{k}>\mathcal{E}^{\mathrm{c}}_{k-q}+\mathcal{E}_{q}^{\mathrm{c}}$, it is possible to find a superposition of excitations with total momentum $k$ whose total energy is less than $\mathcal{E}_{k}$. To get the lower bound on the excitation continuum, we need to compute

[TABLE]

For the system in one dimension, it is sufficient to compute the slope for $k=0$, namely $\mathcal{E}^{\prime}_{0}$, to find the explicit form of the lower bound to be given by

[TABLE]

We can derive the condition on $\mu$, such that there is some range of $k$, for which $\mathcal{E}_{k}^{\min}$ lies underneath the Bogoliubov spectrum $\mathcal{E}^{\mathrm{c}}_{k}$. This condition is given by $\partial^{2}_{k}\mathcal{E}^{\mathrm{c}}_{k}(\mu)=0$, i.e., the second derivative of $\mathcal{E}_{k}^{\mathrm{c}}$ must vanish. Its solution is

[TABLE]

which only exists for $\mu\leq 4$. In higher dimensions, we can consider the slice $k=(k_{x},0,\dots,0)$, which leads to an effective rescaling of $\mu\to\mu+2(\dim-1)$. In this case, we therefore have the condition $\mu<6-2\dim$ to have part of the continuum spectrum to lie underneath the 1-particle dispersion relation $\mathcal{E}_{k}$.

Let us make the following three important remarks. First, the Bogoliubov energy

[TABLE]

is not variational, i.e., it is not the expectation value of the state $\ket{\beta_{0}^{\mathrm{c}},\lambda^{\mathrm{c}}}$ with respect to the full Hamiltonian, but rather the minimal energy of the mean field Hamiltonian $[\hat{H}]_{\ket{\beta^{\mathrm{c}}_{0}}}$. Only the energy $E_{\ket{\beta^{\mathrm{c}}_{0}}}$ is variational, i.e., it minimizes the energy expectation value within the class of coherent states.

Second, the state $\ket{\beta_{0}^{\mathrm{c}},\lambda^{\mathrm{c}}}$ is actually ill defined in the zero mode, due to $\lambda_{k}^{\mathrm{c}}\to\infty$ for $k\to 0$. Put differently, the minimal energy $E_{\mathrm{Bogoliubov}}$ is only reached in the limit of an infinitely squeezed state, whose energy with respect to the full Hamiltonian actually diverges.

Third, we could have computed the Bogoliubov dispersion relation without defining the mean field Hamiltonian $[\hat{H}]_{\ket{\beta^{\mathrm{c}}_{0}}}$, but rather just by studying the real time flow of the full Hamiltonian projected on the manifold of coherent states and linearized around the stationary coherent state $\ket{\beta^{\mathrm{c}}_{0}}$ (see Appendix B). Thus, we do not need to perform the second and third step, if we are only interested in the variational ground state energy $E_{\ket{\beta^{\mathrm{c}}_{0}}}$ and the dispersion relation $\mathcal{E}^{\mathrm{c}}_{k}$, i.e., we are content with not computing $\Delta$. In this case, it is sufficient to linearize Hamiltonian flow around the variational state $\ket{\beta^{\mathrm{c}}_{0}}$ and compute the spectrum of its generator. This is exactly what we do in this paper, but for the extended variational family of all Gaussian states.

Appendix B Bogoliubov theory as coherent TDVP

In this section, we show explicitly that coherent TDVP gives rise to the same 1-particle spectrum as Bogoliubov theory. We consider the manifold of displaced vacua (coherent states)

[TABLE]

Here, $\beta$ is a vector written in the momentum basis, i.e., its components are labeled by $k$. The tangent plane at the state $\ket{\beta}$ is spanned by vectors of the form $\mathcal{D}(\beta)\hat{b}_{k}^{\dagger}\ket{0}$. Therefore, the projected real time evolution can be computed from the quantity

[TABLE]

which for the Bose Hubbard model evaluates to

[TABLE]

Expressed in real components, the resulting evolution is

[TABLE]

where we introduced the transformation matrix

[TABLE]

The symbols $\beta$ and $h$ denote the column vectors that group the values of $\beta_{k}$ and $h_{k}$ for all values of $k$.

The linearization around $\ket{\beta^{\mathrm{c}}_{0}}$ is then given by

[TABLE]

where we expressed the derivatives with respect to $\mathrm{Re}\beta$ and $\mathrm{Im}\beta$ in terms of derivatives with respect to $\beta$ and $\beta^{*}$, taken as independent variables. All derivatives are evaluated at $\beta_{k}=\delta_{k,0}\beta^{\mathrm{c}}_{0}$. The matrix $\mathrm{i}K$, whose eigenvalues $\pm\omega$ represent the TDVP estimate of the 1-particle excitation energies of the model, is then, up to similarity transformations, equal to

[TABLE]

This matrix decomposes into blocks of the form

[TABLE]

with $U^{\mathrm{c}}_{k}=\varepsilon_{k}-2\varepsilon_{0}$ and $V^{\mathrm{c}}_{k}=-\varepsilon_{0}$ as in Appendix A. The eigenvalues $\pm\omega_{k}$ are given by

[TABLE]

which is in full agreement with Bogoliubov theory (34).

Appendix C Iterated Bogoliubov theory

Bogoliubov theory is not self-consistent, in the sense that we construct the mean field Hamiltonian $[\hat{H}]_{\ket{\beta^{\mathrm{c}}_{0}}}$ from the displaced state $\ket{\beta^{\mathrm{c}}_{0}}$, but when we compute the ground state of $[\hat{H}]_{\ket{\beta^{\mathrm{c}}_{0}}}$, we do not find $\ket{\beta^{\mathrm{c}}_{0}}$ again. We can therefore ask if there exists a state $\ket{\psi}$, which is the ground state of its own mean field Hamiltonian $[\hat{H}]_{\ket{\psi}}$.

A natural way to find this state consists of applying the Bogoliubov procedure repeatedly:

Step 1 (displacement). Starting with the Gaussian state $\ket{\beta_{0},\lambda}$ with real $\beta_{0}$ and $\lambda$, we can choose a new real displacement $\beta_{0}^{\prime}$, such that the energy is minimized. The resulting $\beta_{0}^{\prime}$ is given by

[TABLE]

where $u_{k}=\cosh{\lambda}$ and $v_{k}=\sinh{\lambda}$.

Step 2 (quadratic expansion). Once the new displaced ground state $\ket{\beta_{0}^{\prime},\lambda}$ has been found, we can compute the unique quadratic Hamiltonian $[\hat{H}]_{\ket{\beta_{0}^{\prime},\lambda}}$ through the following steps: First, we express the Hamiltonian in terms of the new annihilation operators

[TABLE]

that annihilate $\ket{\beta^{\prime}_{0},\lambda}$. Second, we use the canonical commutation relations $[\delta\hat{b}^{\prime}_{k},\delta\hat{b}^{\prime}_{p}{}^{\dagger}]=\delta_{k,p}$ to write the Hamiltonian as sum of normal ordered operators. Third and finally, we truncate the resulting Hamiltonian at quadratic order to define

[TABLE]

Note that this is the straight-forward generalization of how we constructed the mean field Hamiltonian $[\hat{H}]_{\ket{\beta_{0}^{\mathrm{c}}}}$ in Bogoliubov theory. The coefficients $\tilde{U}_{k}$ and $\tilde{V}_{k}$ are explicitly given by

[TABLE]

Step 3 (new ground state). The ground state $\ket{\beta_{0}^{\prime},\lambda^{\prime}}$ of the mean field Hamiltonian can be encoded in another Bogoliubov transformation

[TABLE]

Here, we have $u_{k}^{\prime}=u_{k}\,\tilde{u}^{\prime}_{k}+v_{k}\,\tilde{v}_{k}$ and $v_{k}^{\prime}=u_{k}\tilde{v}_{k}+v_{k}\tilde{u}^{\prime}_{k}$ with $\tilde{u}^{\prime}_{k}=\cosh{\tilde{\lambda}^{\prime}_{k}}$, $\tilde{v}_{k}=\sinh{\tilde{\lambda}^{\prime}_{k}}$ and $\tanh{2\tilde{\lambda}^{\prime}_{k}}=-\tilde{V}_{k}/\tilde{U}_{k}$. The ground state energy of the mean field Hamiltonian is given by

[TABLE]

with $\tilde{\Delta}^{\prime}=-\sum_{k}\left(\tilde{U}_{k}\,(\tilde{v}_{k}^{\prime})^{2}+\tilde{V}_{k}\,\tilde{u}_{k}^{\prime}\tilde{v}_{k}^{\prime}\right)$. The resulting excitation spectrum is given by

[TABLE]

We can repeat these steps to move from $\ket{\beta_{0}^{\prime},\lambda^{\prime}}$ to $\ket{\beta_{0}^{\prime\prime},\lambda^{\prime}}$, $\ket{\beta_{0}^{\prime\prime},\lambda^{\prime\prime}}$ and so on. For the Bose-Hubbard model, this algorithm converges to the best Gaussian state

[TABLE]

i.e., iterated Bogoliubov theory gives the same approximate ground state as imaginary time evolution on the variational class of all Gaussian states. In fact, we can run into similar troubles as for imaginary time evolution: We could get stuck in a local minima, i.e., a state that is the ground state of its mean field Hamiltonian without being the global minimum among all Gaussian states. Another problem of iterated Bogoliubov theory could arise if we encounter $V^{(n)}_{k}>2U^{(n)}_{k}$ at some step $n$, in which case the mean field Hamiltonian $[\hat{H}]_{\ket{\psi^{(n)}}}$ would not be bounded from below and consequently we could not compute the next ground state.

At $\ket{\psi_{\mathrm{g}}}$, we find the quadratic Hamiltonian

[TABLE]

with $\delta\hat{B}_{k}^{\mathrm{g}}=\lim_{n\to\infty}\delta\hat{B}_{k}^{(n)}$, that provides a self-consistent generalization of traditional Bogoliubov theory. In particular, we find $\ket{\psi_{\mathrm{g}}}$ is the ground state of its own quadratic Hamiltonian $[\hat{H}]_{\ket{\psi_{\mathrm{g}}}}$. Moreover, the ground state energy of $[\hat{H}]_{\ket{\psi_{\mathrm{g}}}}$ coincide with the expectation value $\braket{\psi_{\mathrm{g}}}{\hat{H}}{\psi_{\mathrm{g}}}$ or, put differently, we find $\lim_{n\to\infty}\Delta^{(n)}=0$. However, when looking at the spectrum $\mathcal{E}_{k}^{\mathrm{g}}=\lim_{n\to\infty}\mathcal{E}_{k}^{(n)}$ around the Gaussian ground state approximation $\ket{\psi_{\mathrm{g}}}$, we find a gap $\mathcal{E}_{0}^{\mathrm{g}}>0$, i.e., the 1-particle spectrum computed from iterated Bogoliubov theory does not capture the massless Goldstone mode and is therefore worse than traditional Bogoliubov theory. However, we can use the quadratic Hamiltonian $[\hat{H}]_{\ket{\psi_{\mathrm{g}}}}$ to compute a 2-particle spectrum based on the assumption that the particles with dispersion relation $\mathcal{E}_{k}^{\mathrm{g}}$ do not interact. When comparing the resulting 2-particle continuum with the one from Gaussian TDVP, we find good agreement–in particular in the region around $k=0$.

Appendix D Computation of Gaussian ground state approximation

We review the underlying analytical and semi-analytical methods associated to Section II, that enabled us to compute the best Gaussian state, i.e., the Gaussian state $\ket{\psi_{\mathrm{g}}}$ with the lowest energy expectation value with respect to $\hat{H}$. We can restrict ourselves to searching for the ground state in the translationally invariant submanifold, which is parametrized by $\beta_{0}$ and $\lambda_{k}:=\lambda_{0,k}$. Due to the $\mathrm{U}(1)$ invariance, it is always possible to find a ground state in which both these parameters are real. It turns out to be very convenient to parametrize the state in terms of the Bogoliubov parameters $\beta_{k}$, $u_{k}=\cosh\lambda_{k}$ and $v_{k}=\sinh{\lambda_{k}}$ such that

[TABLE]

will annihilate the Gaussian state $\ket{\beta,\lambda}=\mathcal{U}(\beta,\lambda)\ket{0}$.

The stationary point $\ket{\psi_{\mathrm{g}}}$ of the imaginary time evolution is characterized by vanishing $\mathbb{P}_{\ket{\psi_{\mathrm{g}}}}(-\hat{H})\ket{\psi_{\mathrm{g}}}$, which translates into the conditions

[TABLE]

Rewriting these conditions in terms of our parameters $(\beta_{0},\lambda_{k})$ gives

[TABLE]

where we defined $A=\sum_{k}u_{k}v_{k}$ and $B=\sum_{k}v_{k}^{2}$ and $\varepsilon_{k}$ is the dispersion relation in (23). Equations (84) are solved by

[TABLE]

where we introduced the convenient parameter

[TABLE]

This expression for the solution depends on the final values of the quantities $A$ and $B$ which have to be obtained from the coupled equations

[TABLE]

which can be solved numerically efficiently independently of the dimensionality of the system and with a linear dependence on the system size.

We can similarly express the energy $E$ and particle density $n$ of a Gaussian state $\ket{\psi}$ in terms of $\beta_{0}$, $u_{k}$, $v_{k}$, $A=\sum_{k}u_{k}v_{k}$ and $B=\sum_{k}v_{k}^{2}$ as

[TABLE]

In particular, we can use $\beta_{0}^{\mathrm{g}}$, $u_{k}^{\mathrm{g}}$ and $v_{k}^{\mathrm{g}}$ to compute $E_{\ket{\psi_{\mathrm{g}}}}$.

Appendix E Linearized equations of motion

For calculating the matrix $K$ it is convenient to consider the tangent plane as a real vector space of twice the dimensions compared to the complex one, i.e., we have

[TABLE]

with real linear combinations of $N(N+3)$ tangent vectors. The generator of linearized time evolution is then

[TABLE]

where the matrix $T$ defined in (49), taking the subdivision into blocks to refer to the split between real and imaginary vectors of (92).

The matrix $S$ is block diagonal with each block $S_{k}$ referring to a fixed total momentum. Each block can be written as the sum of a diagonal matrix and a rank 5 matrix, that is $S_{k}=E+CR$, with

[TABLE]

with newly introduced parameters

[TABLE]

where $u_{k}$ and $v_{k}$ have to be evaluated at the solutions corresponding to the ground state approximation defined in Appendix D, i.e., at $\beta_{0}^{\mathrm{g}}$, $u^{\mathrm{g}}_{k}$ and $v^{\mathrm{g}}_{k}$ from (85-87), and $\varepsilon_{k}$ is the dispersion relation in (23).

Given the simple structure of the blocks $S_{k}$, it is easy to diagonalize them numerically. Their eigenvalues are the zeros of the function $f(\omega):=\det[1+R{(E-\omega)}^{-1}C]$. Evaluation only scales linearly with the system size $N$ and moreover, we can characterize analytically some properties of the spectrum in the thermodynamic limit.

More specifically, the function $f(\omega)$ presents a series of poles, given by the diagonal elements of $E$. Its zeros (i.e., the eigenvalues of the system) are positioned one in between each pair of subsequent poles. One subset of the poles, that is the diagonal elements of $\Delta$, for $N\to\infty$ come closer together, creating in the thermodynamic limit a continuous line. The zeros that are in between such poles will therefore also come together to a continuum that represents the continuum in the spectrum of $S_{k}$. The boundaries of this continuum can thus be inferred by computing the values of the minimal and maximal diagonal elements of $\Delta$. In particular, we can identify the minimum, given by $2E_{k}$, with the Higgs excitation mode. In order to give an expression for the Higgs gap at zero momentum, we need to evaluate $2E_{0}$.

We want to do this at constant filling $n$, which is equivalent to imposing $\varepsilon_{0}=-U\left(n+\frac{A+B}{N}\right)$, due to (91). Substituting this condition into equation (90), we find equations for $A$ and $B$ at fixed $n$. These equations admit constant solutions in the limit $U\to 0$. Inserting these solutions in the expression for $E_{0}$, we find the asymptotics of the Higgs gap at constant density for $U\to 0$

[TABLE]

The function $\alpha$ has a complicated analytical expression that admits the large $N$ asymptotics

[TABLE]

Appendix F Linear response theory

We are interested in computing the linear variation $\delta V(t)=\frac{d}{d\lambda}\bra{\psi_{\lambda}(t)}\hat{V}\ket{\psi_{\lambda}(t)}_{\lambda=0}$ due to the perturbation $\lambda\varphi(t)\hat{V}$ of the Hamiltonian. With respect to a set of real coordinates $x^{a}$ of our variational manifold and the corresponding basis $\ket{V_{a}}=\frac{\partial}{\partial x^{a}}\ket{\psi(x^{\mathrm{g}})}$ (where $x^{\mathrm{g}}$ are the coordinates of the stationary point $\ket{\psi_{\mathrm{g}}}$), e.g., the ones of (5) and (92), we find

[TABLE]

with $dV_{a}=\frac{\partial}{\partial x^{a}}\bra{\psi(x)}\hat{V}\ket{\psi(x)}$ and

[TABLE]

Put differently, $\delta\psi^{a}(t)$ is the component of the tangent vector $\ket{\delta\psi(t)}$ in the direction $\ket{V_{a}}$. Note that in (96) as well as in the rest of this appendix we use the Einstein convention for indices, whereby it is understood that all repeated indices are summed over.

Furthermore, we can introduce the components $v^{a}$ of the tangent vector

[TABLE]

representing the linear perturbation of $\ket{\psi_{\mathrm{g}}}$ due to $\hat{V}$.

The key result of the following paragraph is

[TABLE]

where $d\Phi_{t}$ is the linearized flow around the stationary point $\ket{\psi_{\mathrm{g}}}$, often also referred to as push-forward map. The linear response $\delta\psi^{a}(t)$ at time $t$ can be understood as functional $\delta\psi^{a}(t)[\varphi]$ of $\varphi:[-\infty,t]\to\mathbb{R}$ that modulates the perturbation $\hat{V}$ for previous times. We can therefore write

[TABLE]

which states that the $\delta\psi^{a}(t)$ is the superposition of all linear responses

[TABLE]

due to a perturbation with $\varphi(t)=\delta(t-t^{\prime})$. This scenario can be evaluated explicitly to be given by

[TABLE]

Here, the state $\ket{\psi_{\mathrm{g}}}$ is unaffected until time $t^{\prime}$. At this point, it instantaneously kicked to the new state $\ket{\psi_{\lambda}(t^{\prime})}=\Phi^{\hat{V}}_{\lambda}\ket{\psi_{\mathrm{g}}}=\lambda v^{a}\ket{V_{a}}$, where $\Phi^{\hat{V}}_{\lambda}$ represents the projected time evolution with respect to the Hamiltonian $\lambda\delta(t-t^{\prime})\hat{V}$. After the kick, time evolution continuous to be unperturbed and therefore given by $\Phi_{t-t^{\prime}}$ until the time $t$ that we are interested in. Evaluating the derivative with respect to $\lambda$ in (101) therefore gives

[TABLE]

where $d\Phi_{t-t^{\prime}}$ is the linearization of $\Phi_{t-t^{\prime}}$ and $v^{b}\ket{V_{b}}=\frac{d}{d\lambda}|_{\lambda=0}\Phi^{\hat{V}}_{\lambda}\ket{\psi_{\mathrm{g}}}$. Plugging (103) into (100) gives (99).

By observing that

[TABLE]

we finally have the expression

[TABLE]

where we represent the symplectic form as

[TABLE]

Next, we consider the Fourier transformed response function, that can be written as

[TABLE]

We define $Z_{V}(\omega)$ as the spectral response function.

In order to obtain the expression (109), it is useful to give a spectral decomposition of the linear operator $d\Phi_{t}$. As mentioned, if we consider as the initial state of our evolution the approximate ground state $\ket{\psi_{\mathrm{g}}}$, $d\Phi_{t}$ becomes a linear map from the tangent plane $\mathcal{T}_{\ket{\psi_{\mathrm{g}}}}\mathcal{M}$ onto itself, given by the exponential of the generator $K$, i.e., we have

[TABLE]

The generator $K$ was computed in Appendix E and, as discussed before, its eigenvectors $e^{a}(\omega)$ appear in complex conjugate pairs satisfying

[TABLE]

This means that $K$ can be decomposed as

[TABLE]

with

[TABLE]

where the $\omega_{i}$ are all taken to be positive.

Using this, the object in the square brackets in (105) can be written as

[TABLE]

where we introduced the terms

[TABLE]

such that the vectors $e(\omega)$ are normalized to satisfy $\delta(\omega)=\pm 1$.

Taking the Fourier transform as in (109), we find

[TABLE]

The generator $K^{a}{}_{b}$ is a linear map on the tangent space, taken as real vector space. It is not complex-linear, i.e., $K^{a}{}_{b}$ does not commute with $J^{a}{}_{b}$, which is the linear map representing multiplication with the imaginary unit

[TABLE]

In particular, this implies that $K$ is not anti-Hermitian with respect to the standard inner product, which means that $d\Phi_{t}=e^{tK}$ is not unitary.

However, it is easy to find an alternative inner product on the tangent space that turns $K$ into an anti-Hermitian operator and is therefore preserved by the linearized flow $d\Phi_{t}$. Given two vectors $\ket{x}=x^{a}\ket{V_{a}}$ and $\ket{y}=y^{a}\ket{V_{a}}$, we can define the inner product

[TABLE]

which is a well-defined Hermitian inner product. With respect to this inner product, $K$ is anti-Hermitian, i.e., $K^{\dagger}=-K$, which implies that $d\Phi_{t}=e^{tK}$ is unitary.

The modified inner product (124) on the tangent plane can be used to calculate overlaps between vectors, consistent with the unitary time evolution. In particular, the ovelaps between the evolving perturbation and the 1-particle subspace of the tangent plane shown in Figure 8 are calculated with this inner product.

Appendix G Random phase approximation

In this section, we show the fluctuation spectrum from the Gaussian state approach can also be obtained from the standard random phase approximation (RPA) based on the ladder diagram, similarly to what has been done for fermions in Demler et al. (1996).

Around any Gaussian variational state $\ket{\psi}$, characterized by the real numbers $\beta_{0}$, $u_{k}$ and $v_{k}$ as described in Appendix D, the Hamiltonian $H=E+\sum_{j=1}^{4}h_{j}$ is decomposed into five normal ordered terms (with respect to the given Gaussian state $\ket{\psi}$) via the Wick theorem. We will here be interested in taking this reference state to be our self-consistent Gaussian ground state approximation $\ket{\psi_{\mathrm{g}}}$ parametrized by real $\beta_{0}^{\mathrm{g}}$, $u_{k}^{\mathrm{g}}$ and $v^{\mathrm{g}}_{k}$ from (85-87).

For this particular state, the energy $E$ is given by (91). The linear term $h_{1}$ vanishes. The quadratic term, when written in terms of the Bogoliubov-transformed operators

[TABLE]

takes the diagonal form

[TABLE]

where the dispersion relation $\mathcal{E}_{k}=\sqrt{U_{k}^{2}-V_{k}^{2}}$ with

[TABLE]

coincides with (77).

The cubic and quartic terms become

[TABLE]

We calculate the time-ordered correlation functions

[TABLE]

by treating $h_{3}$ and $h_{4}$ as perturbations. To the 1-loop order, as shown in Figure 12, the correlation functions read

[TABLE]

where the free propagators for single and two excitations are $G_{0}(k,\omega)=(\omega-\mathcal{E}_{k}+\mathrm{i}0^{+})^{-1}$ and $\Pi_{0}(k,p,\omega)=(\omega-\mathcal{E}_{k+p}-\mathcal{E}_{p}+\mathrm{i}0^{+})^{-1}$.

The coupling between single and two excitations is described by interaction vertices

[TABLE]

and the interaction vertices of two excitatations are

[TABLE]

In the compact matrix form, we can read (134-137) as

[TABLE]

in the basis $\mathbf{G}=(G(k),F(k),\Pi(k,p),\bar{\Pi}(k,p))^{\intercal}$, where $\mathbf{V}=(1,0,0,0)^{\intercal}$, and the matrix

[TABLE]

It can be verified that, if we evaluate the parameters $\beta_{0}$, $u_{k}$ and $v_{k}$ to be the ones corresponding to the best Gaussian state $\ket{\psi_{\mathrm{g}}}$, the matrix $\mathbf{M}$ is exactly the linearized time evolution generator matrix $K$ transfromed into the Bogoliubov basis, whose eigenvalues are the poles of the Green function $\mathbf{G}$ thus determine the fluctuation spectrum.