Flow equations for cold Bose gases

TL;DR

This paper develops flow equations for cold Bose gases to accurately compute ground state energies, using a self-consistent truncation scheme and benchmarking against the Lieb-Liniger model.

Contribution

It introduces a novel flow equation approach with a self-consistent truncation for three-body operators in cold Bose gases, enabling precise energy calculations.

Findings

Good agreement with Lieb-Liniger model for weak/moderate interactions

Effective estimation of truncation errors

Method applicable to general N-body systems

Abstract

We derive flow equations for cold atomic gases with one macroscopically populated energy level. The generator is chosen such that the ground state decouples from all other states in the system as the renormalization group flow progresses. We propose a self-consistent truncation scheme for the flow equations at the level of three-body operators and show how they can be used to calculate the ground state energy of a general $N$-body system. Moreover, we provide a general method to estimate the truncation error in the calculated energies. Finally, we test our scheme by benchmarking to the exactly solvable Lieb-Liniger model and find good agreement for weak and moderate interaction strengths.

Equations

$$\mathbf{M}(\mathbf{Q})=\mathbf{M}+[\mathbf{\eta },\mathbf{M}]\delta s+\dots .$$

$$\frac{\mathrm{d}{M}_{ij}}{\mathrm{d}s}=\sum _k({\eta }_{ik}{M}_{kj}-M_{ik}{\eta }_{kj}).$$

$$\frac{\mathrm{d}{M}_{ij}}{\mathrm{d}s}=\alpha (M_{ib}{\delta }_{ja}+M_{jb}{\delta }_{ia}-M_{ia}{\delta }_{bj}-M_{ja}{\delta }_{bi}),$$

$$\frac{\mathrm{d}{M}_{ab}}{\mathrm{d}s}=\frac{\mathrm{d}{M}_{ba}}{\mathrm{d}s}=\alpha (M_{bb}-M_{aa}).$$

$$\frac{\mathrm{d}{M}_{0i>0}}{\mathrm{d}s}=-M_{00}{\alpha }_i+\sum _k{\alpha }_k{M}_{ki},$$

$$\frac{\mathrm{d}{M}_{i>0j>0}}{\mathrm{d}s}=-M_{0j}{\alpha }_i-M_{i0}{\alpha }_j.$$

$$\frac{\mathrm{d}{M}_{0i>0}}{\mathrm{d}s}=-\sum _{k\ne 0}(M_{ik}-M_{00}{I}_{ik})M_{0k}.$$

$$M_{0i}(s)=\sum _{k\ne 0}{\left(\text{mathbfcal}T{e}^{-{\int }_0^s(\mathbf{M}-M_{00}\mathbf{I}{)}^{\prime }\mathrm{d}s}\right)}_{ik}{M}_{0k}(0),$$

$$\mathbf{H}=A_{ij}{\mathbf{a}}_i^{†}{\mathbf{a}}_j+\frac{1}{2}{B}_{ijkl}{\mathbf{a}}_i^{†}{\mathbf{a}}_j^{†}{\mathbf{a}}_k{\mathbf{a}}_l,$$

$$:{\mathbf{a}}_{{\alpha }_1}^{†}{\mathbf{a}}_{{\alpha }_2}:$$

$$:{\mathbf{a}}_{{\alpha }_1}^{†}{\mathbf{a}}_{{\alpha }_2}^{†}{\mathbf{a}}_{{\alpha }_3}{\mathbf{a}}_{{\alpha }_4}:$$

$$\mathbf{H}=ϵN\mathbf{I}+f_{ij}:{\mathbf{a}}_i^{†}{\mathbf{a}}_j:+\frac{1}{2}{\Gamma }_{ijkl}:{\mathbf{a}}_i^{†}{\mathbf{a}}_j^{†}{\mathbf{a}}_k{\mathbf{a}}_l:,$$

$$ϵN$$

$$f_{{\alpha }_1{\alpha }_2}$$

$${\Gamma }_{{\alpha }_1{\alpha }_2{\alpha }_3{\alpha }_4}$$

$$\frac{\mathrm{d}\mathbf{H}(s)}{\mathrm{d}s}=[\mathbf{\eta }(s),\mathbf{H}(s)],$$

$$\mathbf{\eta }(s)$$

$$\mathbf{H}(s)$$

$$:{\mathbf{a}}_{{\alpha }_1}^{†}{\mathbf{a}}_{{\alpha }_2}^{†}{\mathbf{a}}_{{\alpha }_3}{\mathbf{a}}_{{\alpha }_4}:=:{\mathbf{a}}_{{\alpha }_2}^{†}{\mathbf{a}}_{{\alpha }_1}^{†}{\mathbf{a}}_{{\alpha }_3}{\mathbf{a}}_{{\alpha }_4}:=:{\mathbf{a}}_{{\alpha }_1}^{†}{\mathbf{a}}_{{\alpha }_2}^{†}{\mathbf{a}}_{{\alpha }_4}{\mathbf{a}}_{{\alpha }_3}:.$$

$$[\mathbf{\eta },\mathbf{H}{]}^{(3)}=({\eta }_{iklj}{\Gamma }_{jbcd}-{\Gamma }_{iklj}{\eta }_{jbcd}){\mathbf{a}}_i^{†}{\mathbf{a}}_k^{†}{\mathbf{a}}_b^{†}{\mathbf{a}}_l{\mathbf{a}}_c{\mathbf{a}}_d,$$

$$\frac{\mathrm{d}ϵ}{\mathrm{d}s}=S_{00}+(N-1)\left(\frac{1}{2}{S}_{00ii00}-S_{000000}\right),$$

$$\frac{\mathrm{d}{f}_{{\alpha }_1{\alpha }_2}(s)}{\mathrm{d}s}=-(N-1{)}^2(S_{0{\alpha }_{1}00{\alpha }_{2}0}+S_{0{\alpha }_1{\alpha }_{2}000}+S_{000{\alpha }_1{\alpha }_{2}0})+(N-1)S_{{\alpha }_{1}0ii0{\alpha }_2}$$

$$+S_{0{\alpha }_1{\alpha }_{2}0}+\frac{(N-1)(N-2)}{2}(S_{00{\alpha }_2{\alpha }_{1}00}+S_{0{\alpha }_{1}00{\alpha }_{2}0}{D}_{{\alpha }_2}{D}_{{\alpha }_1})$$

$$+\frac{(N-1)(N-2)}{2}(S_{0{\alpha }_1{\alpha }_{2}000}{D}_{{\alpha }_1}+S_{000{\alpha }_1{\alpha }_{2}0}{D}_{{\alpha }_2})+S_{{\alpha }_1{\alpha }_2},$$

$$\displaystyle\frac{\mathrm{d}\Gamma_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}(s)}{\mathrm{d}s}=\frac{(1+P_{\alpha_{1}\alpha_{2}})(1+P_{\alpha_{3}\alpha_{4}})}{2}\bigg{(}S_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}$$

$$-(N-1)(S_{{\alpha }_1{\alpha }_2{\alpha }_{3}00{\alpha }_4}+S_{{\alpha }_{1}00{\alpha }_2{\alpha }_3{\alpha }_4}+\frac{1}{2}{S}_{{\alpha }_1{\alpha }_{2}ii{\alpha }_3{\alpha }_4})$$

$$+(N-2)D_{{\alpha }_1}{D}_{{\alpha }_4}{S}_{0{\alpha }_1{\alpha }_3{\alpha }_2{\alpha }_{4}0}+(N-2)I_{{\alpha }_1{\alpha }_2}{D}_{{\alpha }_4}{S}_{{\alpha }_1{\alpha }_2{\alpha }_{3}00{\alpha }_4}$$

$$\displaystyle+(N-2)D_{\alpha_{1}}I_{\alpha_{3}\alpha_{4}}S_{\alpha_{1}00\alpha_{2}\alpha_{3}\alpha_{4}}+(N-2)I_{\alpha_{1}\alpha_{2}}I_{\alpha_{3}\alpha_{4}}S_{\alpha_{1}\alpha_{2}ii\alpha_{3}\alpha_{4}}\bigg{)},$$

$$S_{{\alpha }_1{\alpha }_2}^{(1)}$$

$$S_{{\alpha }_1{\alpha }_2{\alpha }_3{\alpha }_4}^{(2)}$$

$$S_{{\alpha }_1{\alpha }_2{\alpha }_3{\alpha }_4{\alpha }_5{\alpha }_6}$$

$$\frac{\mathrm{d}{\mathbf{H}}_a(s)}{\mathrm{d}s}$$

$${\mathbf{H}}_a(s)={\int }_0^s\mathbf{W}\mathrm{d}x+{\int }_0^s[\mathbf{\eta }(x),{\mathbf{H}}_a(x)]\mathrm{d}x,$$

$${\mathbf{H}}_a(s)=\mathbf{U}(s)\left({\int }_0^s{\mathbf{U}}^{†}(x)\mathbf{W}\mathbf{U}\mathrm{d}x\right){\mathbf{U}}^{†}(s),$$

$$\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}s}=\mathbf{\eta }\mathbf{U}\to \mathbf{U}(s)=\text{mathbfcal}T{e}^{{\int }_0^s\mathbf{\eta }(x)\mathrm{d}x}.$$

$$-\frac{1}{2}\sum _{\alpha =1}^N\frac{{\partial }^2}{\partial x_{\alpha }^2}\Psi +g\sum _{{\alpha }_1<{\alpha }_2}\delta (x_{{\alpha }_1}-x_{{\alpha }_2})\Psi =E_N\Psi ,$$

$$\mathbf{\eta }(s)=f_{i0}(s):{\mathbf{a}}_i^{†}{\mathbf{a}}_0:+\frac{1}{2}{\Gamma }_{ij00}(s):{\mathbf{a}}_i^{†}{\mathbf{a}}_j^{†}{\mathbf{a}}_0{\mathbf{a}}_0:-H.c..$$

$$\delta e\simeq \frac{1}{N}\sum _p\frac{{\left(⟨{\Phi }_p|{\int }_0^{\infty }\mathbf{W}(s)\mathrm{d}s|\Phi ⟩\right)}^2}{⟨\Phi |\mathbf{H}|\Phi ⟩-⟨{\Phi }_p|\mathbf{H}|{\Phi }_p⟩},$$

Full text

Keywords: cold Bose gases, similarity renormalization group, mesoscopic systems

1 Introduction

The worlds of many- and few-body physics are generally far apart. In the former, the number of particles is often infinite, while few-body systems normally do not contain more than a handful of particles. The typical goal in many-body physics is to calculate thermodynamic quantities such as the energy per particle and the density profile. However, the large number of degrees of freedom in many-body systems usually means that various approximations and/or large computational resources are needed to achieve this goal. In contrast, it is often possible to solve few-body problems exactly, i.e., to find the full spectrum of the Hamiltonian and the corresponding wave functions. There is an interesting class of systems that are in between these two extremes. These are finite systems in which the number of particles is sufficiently large for many-body phenomena, such as superfluidity or Bose-Einstein condensation, to emerge [1, 2, 3]; but they are still small enough to be within reach for numerically exact ab-initio calculations that use microscopic Hamiltonians. The investigation of such finite systems is crucial to understand how many-body phenomena arise from few-body body physics and microscopic interactions of the constituents.

To investigate this progression from few- to many-body behavior theoretically one needs reliable numerical techniques in the transition region. There exists a number of suitable techniques in physics and chemistry and new methods are being developed (see, e.g., references [4, 5, 6, 7, 8, 9, 10]). A significant breakthrough was made with the development of flow equation methods which are also referred to as the similarity renormalization group (SRG) [11, 12]. In this approach, a set of differential equations is solved to obtain unitarily equivalent Hamiltonians with desirable properties. This set is determined by a generator, which controls the change of the Hamiltonian at every step of the evolution. Note that this generator is determined dynamically. Its matrix elements depend on the flow parameter $s$ and are calculated at every step of the evolution from the transformed Hamiltonian. This represents one of the key advantages of the SRG, which allows one to find a (block-)diagonal representation of the Hamiltonian.

Recently a new approach based on flow equations, the in-medium similarity renormalization group (IMSRG), has been proposed for nuclear physics problems where the fundamental degrees of freedom are fermionic [13]. The IMSRG is a very promising method for medium-mass nuclei, which lie exactly in the few- to many-body transition region discussed above (see [14] for a recent review).

In this paper, we develop a similar method for cold Bose gases. To this end, we write flow equations for bosonic systems with a macroscopic occupation of one state. We introduce a suitable truncation scheme that facilitates numerical calculations and discuss its accuracy. In particular, we provide an algorithm to estimate the truncation error using perturbation theory. We validate our method using the exactly solvable Lieb-Liniger model in one dimension, and show that even without preliminary knowledge of the reference state our method can be used to accurately describe systems with weak and intermediate interaction strengths.

The paper is organized as follows: in section 2, we review the foundations of the SRG method. In section 3, we introduce the Hamiltonian of interest and write down the flow equations to find its eigenvalue. Here we also discuss the accuracy of our approach and provide a way to estimate the accuracy of the calculated energies. We test the method in section 4 using the exactly solvable Lieb-Liniger model as a benchmark. Section 5 concludes the paper with a summary of our results and a brief outlook on the generalization to three spatial dimensions. For the reader’s convenience, we include six appendices with technical details on the evaluation of commutators, the truncation of the three-body operator, the convergence of the two-body energy, the effective interaction used in the Lieb-Liniger model, the use of White-type generators, and the error estimation.

2 Preliminaries

For a self-contained discussion, we first review the SRG method as it forms the basis of our approach (cf. [11, 12, 14, 15, 16]). To this end, we introduce a real symmetric matrix ${\bf M}$ that represents a linear operator in a particular basis111Two comments are in order here. First, we use bold type for matrices and operators, e.g., ${\bf M}$, and italic type for the corresponding matrix elements, e.g., $M_{ij}$. Second, we choose to work with a real matrix ${\bf M}$ to simplify the discussion. The ideas presented here can be extended straightforwardly to Hermitian matrices.. If we transform this basis using some orthogonal matrix ${\bf Q}$ (i.e., ${\bf QQ}^{T}={\bf I}$, where ${\bf I}$ is the identity matrix) then the linear operator will be represented by the new matrix ${\bf M}({\bf Q})\equiv{\bf QMQ}^{T}$, which is unitarily equivalent to $\mathbf{M}$. The SRG equations simply describe the change of ${\bf M}$ for a small change of the basis: ${\bf Q}={\bf I}+{\boldsymbol{\eta}}\delta s$ ($|\delta s|\ll 1$),

[TABLE]

In the limit $\delta s\to 0$, the SRG equations can be written in the differential form:

[TABLE]

They define the evolution of matrix elements $M_{ij}$ driven by the skew-symmetric matrix ${\boldsymbol{\eta}}=-{\boldsymbol{\eta}}^{T}$. By specifying ${\boldsymbol{\eta}}$, one finds a unitarily equivalent to ${\bf M}$ matrix with some desired properties. Note that the system of equations (2) is often called the “flow equation”, as it defines the “flow” of matrix elements under the SRG transformation, and the “generator” ${\boldsymbol{\eta}}$ determines the flow by defining the “direction” of the transformation at each value of $s$.

We illustrate the evolution using a generator ${\boldsymbol{\eta}}$ that contains only two non-zero elements $\eta_{ab}=-\eta_{ba}$, i.e., $\eta_{ik}(s)=\alpha(s)(\delta_{ia}\delta_{kb}-\delta_{ib}\delta_{ka})$. This matrix leads to the system of equations

[TABLE]

in which the element $M_{ab}=M_{ba}$ is transformed as

[TABLE]

Let us assume that we want the flow to eliminate the element $M_{ab}$ as $s\to\infty$, e.g., by demanding that $M_{ab}(s)=M_{ab}(0)e^{-s}$. Inserting this ansatz into (4) we find that $\alpha=-M_{ab}/(M_{bb}-M_{aa})$ fulfills this requirement222Note that to eliminate $M_{ab}$, we could also have chosen $\alpha(s)=-M_{ab}f(s)/(M_{bb}-M_{aa})$ with $f(s)>0$, e.g., $f(s)=|M_{bb}-M_{aa}|$, as then $\frac{\mathrm{d}M^{2}_{ab}}{\mathrm{d}s}<0$ if $M_{ab}\neq 0$, which means that $M_{ab}(s)$ dies off., i.e., it decouples the basis states with numbers $a$ and $b$. Note, however, that to achieve this decoupling, the flow usually needs to couple states that were not coupled before. For example, if we had $M_{cd}=0$ at $s=0$, then this element will attain a non-zero value if $M_{cb}\delta_{da}+M_{db}\delta_{ca}-M_{ca}\delta_{bd}-M_{da}\delta_{bc}\neq 0$.

Let us give another example of how one can obtain a new matrix with some desired properties by choosing an appropriate generator ${\boldsymbol{\eta}}$. To this end, we use a generator that contains only one row and one column, i.e., $\eta_{ik}=\delta_{i0}\alpha_{k}-\delta_{k0}\alpha_{i}$ with $\alpha_{0}=0$. The corresponding flow equations are

[TABLE]

The prescription $\alpha_{i>0}=-M_{0i}$, which is inspired by the previous example, leads to

[TABLE]

A formal solution to this equation can be found using the Magnus expansion

[TABLE]

here $\mathbfcal{T}$ denotes the $s$-ordering operator (see, e.g., [16, 17]), and ′ means that the first row and the first column should be crossed out from the matrix. If all $M_{0j}$ are initially small (i.e., much smaller than the differences of the eigenvalues of $\bf{M}$), then the long time behavior can be estimated by examining the matrix $({\bf M}(0)-M_{00}(0){\bf I})^{\prime}$. This shows that if $M_{00}(0)$ is close to the ground state then $M_{0i}$ is driven to zero during the evolution, and hence $M_{00}(s\to\infty)$ is the ground state of the matrix. These considerations can be useful in physics problems, as they allow one to find eigenenergies of a system by diagonalizing (block-diagonalizing) the corresponding Hamiltonian. This statement will be exemplified below.

3 Flow Equations

3.1 Hamiltonian

We now consider a system of $N$ bosons that is described by the Hamiltonian

[TABLE]

where ${\boldsymbol{a}_{\alpha_{1}}}$ is the standard annihilation operator333From now on we adopt in the numbered equations the Einstein summation convention for the letters from the Latin alphabet, i.e., $A_{ij}{\boldsymbol{a}}_{i}^{\dagger}{\boldsymbol{a}}_{j}\equiv\sum_{ij}A_{ij}{\boldsymbol{a}}_{i}^{\dagger}{\boldsymbol{a}}_{j}$, and reserve the indices $\alpha_{1,2,...}$ for the places where this convention is not implied.. Since the system is bosonic, ${\bf H}$ is symmetrized with respect to particle exchanges, i.e., $B_{ijkl}=B_{jikl}=B_{ijlk}$. For our numerical calculations this Hamiltonian should be written as a finite-dimensional matrix. Therefore, we assume that the sums in every index run only up to some number $n$ that defines the dimension of the used one body basis.

We are mainly interested in the ground state properties of systems with a macroscopic population of one state (condensate). To incorporate our intentions in the Hamiltonian, we normal order operators using the reference state $\Phi=\prod_{\alpha=1}^{N}f(x_{\alpha})$, where $f(x)$ is some one body function that approximates the condensate (e.g., obtained by solving a suitable Gross-Pitaevski equation):

[TABLE]

where ${\bf P}_{\alpha_{1}\alpha_{2}}$ exchanges the indices $\alpha_{1}$ and $\alpha_{2}$, $\kappa\equiv\frac{N-1}{2N}$, $\rho_{\alpha_{1}\alpha_{2}}\equiv\langle\Phi|{\boldsymbol{a}}_{\alpha_{1}}^{\dagger}{\boldsymbol{a}}_{\alpha_{2}}|\Phi\rangle$, $\rho_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}\equiv\langle\Phi|{\boldsymbol{a}}_{\alpha_{1}}^{\dagger}{\boldsymbol{a}}_{\alpha_{2}}^{\dagger}{\boldsymbol{a}}_{\alpha_{3}}{\boldsymbol{a}}_{\alpha_{4}}|\Phi\rangle$. These operators connect the reference state to the states that contain one and two excitations respectively.

Using the normal-ordered operators we rewrite the Hamiltonian as

[TABLE]

where

[TABLE]

$\epsilon$ is the energy per particle in the reference state, and the elements $f_{\alpha_{1}\alpha_{2}}$ and $\Gamma_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}$ describe one- and two-body excitations, correspondingly. We will construct the Hamiltonian matrix using the basis that contains $f(x)$ as the zero element, therefore, from now on we use $\rho_{\alpha_{1}\alpha_{2}}=\delta_{\alpha_{1}0}\delta_{\alpha_{2}0}N$ and $\rho_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}=\delta_{\alpha_{1}0}\delta_{\alpha_{2}0}\delta_{\alpha_{3}0}\delta_{\alpha_{4}0}N(N-1)$.

3.2 Truncated flow equations

Our goal is to find a matrix representation of ${\bf H}$ in which the couplings to the reference state vanish, i.e., $f_{i0}=\Gamma_{ij00}=0$, so $\epsilon$ is an eigenenergy. To achieve this, we write ${\bf H}$ in a particular basis and then use the flow equations

[TABLE]

where the antihermitian matrix ${\boldsymbol{\eta}}$ eliminates the couplings. To solve this equation, we assume that during the flow the generator and the Hamiltonian contain only one- and two-body operators, i.e.,

[TABLE]

For now we leave the parameters $\xi_{ij}$ and $\eta_{ijkl}$ undetermined. We just mention that they must be chosen such that the couplings vanish at $s\to\infty$. This is usually achieved by calculating $\xi_{ij}(s)$ and $\eta_{ijkl}(s)$ for every $s$ from the evolved matrix elements of the Hamiltonian. We give a possible choice of ${\boldsymbol{\eta}}$ in the next section. It is worthwhile noting that since ${\boldsymbol{\eta}}$ is antihermitian, the following relations must be satisfied $\xi_{ji}=-\xi_{ij}^{*}$, and $\eta_{klij}=-\eta_{ijkl}^{*}$. Moreover, we assume that $\eta_{ijkl}=\eta_{jikl}=\eta_{ijlk}$, because by construction

[TABLE]

Note that equations (16), (17) and (18) do not lead in a general case to a self-consistent system of equations. Indeed, the commutator444From now on we omit the argument $s$ whenever it cannot cause confusion. $[{\boldsymbol{\eta}},{\bf H}]$ contains the three body operator (see A)

[TABLE]

where the superscript $(3)$ corresponds to the piece of the commutator that contains three-body operators. This piece is apparently beyond the scheme put forward in (17) and should be omitted. To this end, we extract from $[{\boldsymbol{\eta}},{\bf H}]^{(3)}$ the terms that contain at least one operator ${\boldsymbol{a}}_{0}^{\dagger}{\boldsymbol{a}}_{0}$, and put to zero the remaining pieces (called ${\bf W}$). The operator ${\boldsymbol{a}}_{0}^{\dagger}{\boldsymbol{a}}_{0}$ is then treated as a constant because of the assumed macroscopic occupation of the lowest state (see B).

After the three-body operator is truncated, we end up with a closed system of equations. To write it down, we equate the coefficients in front of the same operators, i.e.,

[TABLE]

where $D_{\alpha_{1}}=2-\delta_{\alpha_{1}0}$, $I_{\alpha_{1}\alpha_{2}}=1+\delta_{\alpha_{1}0}\delta_{\alpha_{2}0}-2\delta_{\alpha_{2}0}$, and

[TABLE]

This system of equations can be solved using standard solvers of ordinary differential equations. During the evolution, an appropriate choice of ${\boldsymbol{\eta}}$ eliminates the couplings $f_{i0}=0$ and $\Gamma_{ij00}=0$, so that $\epsilon(s\to\infty)$ approximates an eigenenergy of the system. It is worthwhile noting that it is not guaranteed that $\epsilon$ is close to the ground state energy, unless $\Phi$ describes the ground state wave function “well” (so that $f_{i0}$ and $\Gamma_{ij00}$ are much smaller than the differences of the eigenenergies of ${\bf H}$).

3.3 Error estimation

Since the flow equations are truncated at the level of three-body operators and beyond, it is important to estimate the error induced by this approximation. Let us imagine that we have integrated the flow equations (21)-(23) up to $s\to\infty$, and obtained the operator ${\bf H}(s)$ within our truncation scheme as well as the generator ${\boldsymbol{\eta}}(s)$. Now we assume that ${\boldsymbol{\eta}}$ is fixed for every $s$ and use it to introduce the operator $\mathbfcal{H}$ that solves equation (16) with the initial condition $\mathbfcal{H}(s=0)={\bf H}(s=0)$ without any truncations. Hence, $\mathbfcal{H}$ is unitarily equivalent to ${\bf H}(0)$. We emphasize that ${\boldsymbol{\eta}}(s)$ is given from the beginning for every $s$ and not obtained dynamically as before.

The operator $\mathbfcal{H}$ can be written as $\mathbfcal{H}(s)={\bf H}(s)+{\bf H}_{a}(s)$, where ${\bf H}(s)$ is obtained from the truncated flow and ${\bf H}_{a}$ satisfies the equation

[TABLE]

supplemented by the initial condition ${\bf H}_{a}(s=0)=0$. Note that the operator ${\bf H}_{a}$ is generated by ${\bf W}$, which is the part of (20) that is neglected in our truncation scheme. Therefore, we postulate that our approximation is meaningful only if ${\bf H_{a}}(s\to\infty)$ can be treated as a small perturbation for the state of interest. In this case $\epsilon(s\to\infty)$ is close to the exact eigenenergy of the operator $\mathbfcal{H}$.

To estimate ${\bf H_{a}}(s)$ we write two formal solutions to (25)

[TABLE]

where $\bf{U}$ is the transformation matrix generated by ${\boldsymbol{\eta}}$

[TABLE]

Equations (26) and (27) allow us to estimate ${\bf H}_{a}(s\to\infty)$ and then use matrix perturbation theory to find the correction to the energy of the eigenstate. We will illustrate this procedure below using the Lieb-Liniger model.

4 Lieb-Liniger Model

To test our method, we use the exactly solvable Lieb-Liniger model [18], which describes $N$ spinless bosons on a ring of length $L$. The particles interact via delta functions, so the corresponding one-dimensional Schrödinger equation is

[TABLE]

where we put $\hbar=m=1$ for convenience. The parameters of the model are $\gamma=g/\rho$ and $e=2E_{N}/(N\rho^{2})$, where $\rho=N/L$ is the density of the system. Since this model is exactly solvable for any $N,L$ and $g$, it gives us a good reference point for testing our approach. Note, however, that we do not expect our approach to work extremely well for large systems, as strong correlations preclude the existence of a “true” BEC in one spatial dimension.

To write the initial matrix elements and the reference state, we use the one-body basis of plane waves, i.e., $\phi_{i}(x)=e^{ik_{i}x}/\sqrt{L}$, where $k_{i}\in\{0,\pm 1,\pm 2,...\}2\pi/L$ and $\Phi=L^{-N/2}$. Inspired by the discussion in section 2, we write the generator as

[TABLE]

Here we explicitly relate the parameters of the generator (17) to the parameters of the Hamiltonian (18) for every $s$. This generator decouples the element $\langle\Phi|{\bf H}|\Phi\rangle$ from the rest; see figure 1, where we plot $\epsilon(s)$ and $\sum_{p>0}|\langle\Phi_{p}|{\bf H}(s)|\Phi\rangle|^{2}$ with $\Phi_{p}$ containing one- and two-body excitations. Therefore, the latter represents the coupling to the state of interest. We see that during the flow the couplings vanish, and $\epsilon(s\to\infty)$ can be interpreted as the eigenvalue of the matrix. Note that we do not plot any numbers on the $y$ axis as this schematic plot is representative for all considered cases. In our code we use units with $L=2\pi$, which gives a particularly simple form of the momenta, and sets the scale for the energy and $s$ in the problem. For example, the energy difference between the two lowest non-interacting states is one, and therefore the slowest dynamics in the weakly interacting case are described approximately by $e^{-s}$. Figure 1 shows that in these units the decoupling indeed occurs for $s$ of ${\cal O}(1)$ as expected. The study of the flow for other generators is beyond the scope of the present paper. However, we did check (see E) that the results obtained with the generator (30) agree with the results obtained using White’s generator [14, 15], which includes additional energy denominators compared to (30).

4.1 Results

$N=4$. We start with $N=4$. Note that for a cutoff $n\simeq 25$ (in the one-body sector this corresponds to the maximal energy of $288\pi^{2}/L^{2}$), we can easily run the flow until the states are decoupled with high accuracy. Therefore, we have only two sources of error. The first is due to the truncation of the Hamiltonian at $s=0$. This error vanishes in the limit of large $n$, but since the delta function potential has a hard core (it couples all plane waves equally strongly) the convergence to the $n\to\infty$ limit might be relatively slow (see C). However, one can still extract accurate results either by fitting (see C) or by using an effective interaction (see D). To be on the safe side, we first solve the problem using the former method and then using the latter. The results of both methods agree well. This is demonstrated explicitly in figures 5 and 6 for two parameter sets.

The second error is due to the truncation of the three-body term in Eq. (20). To estimate this error, we note that according to (26) for a weak interaction ${\bf H}_{a}\simeq\int\bf{W}\mathrm{d}s$. By definition, the operator $\bf{W}$ connects the state of interest to the states with three excitations. To calculate the contribution to the energy of the perturbation ${\bf H_{a}}$, we use the standard second-order eigenvalue correction from perturbation theory, i.e.,

[TABLE]

where the sum goes over the all states that contain three particles excited out of the condensate. For consistency, we will keep only the lowest terms in $g$ in the denominator.

We show our results in figure 2. On the scale of the figure, the results for the bare delta-function interaction and the effective interaction are indistinguishable. We see that the SRG reproduces the exact results at weak and moderate coupling strengths. However, when the interaction strength increases the energy starts to deviate noticeably. This behavior can be understood by calculating $\delta e$. We see that this term grows very rapidly (numerical analysis reveals that in the considered interval this term grows faster than $\gamma^{2}$) and already at $\gamma=\pi^{2}/2$ it accounts for about $25$% of the SRG result. This shows that the used truncation scheme is not accurate for this $\gamma$ making us stop our calculations.

$N=15$. Our results for $N=15$ are shown in figure 3. On the scale of the figure the results for the bare delta-function interaction and the effective interaction are again indistinguishable. We see a similar trend as for $N=4$: The SRG reproduces well the exact results at small and moderate coupling strength, but fails to describe strongly interacting systems. The window of applicability of the SRG for $N=15$ is slightly smaller than for $N=4$, which is expected from our error estimation which shows that $\delta e$ grows with $N$ (see F).

5 Conclusions

In this paper we have developed a non-perturbative numerical procedure to address bosonic systems with a macroscopic occupation of one state. The method is based on the SRG approach in which the Hamiltonian is transformed to decouple the state of interest from the rest. This transformation is done through a sequence of infinitesimally small rotations in the state space described by a system of differential equations. To make this system solvable with the standard numerical software, we truncate it at the level of three-body operators, and present means to estimate the introduced uncertainty. To illustrate our approach we turn to the Lieb-Liniger model, which shows that our flow equations describe small systems with weak and moderate interactions well. Note that our method can be used to describe two- and three-dimensional systems and we use here a one-dimensional model because its exact solutions allow us to directly test our procedure (although studies of trapped systems in one spatial dimension are interesting on their own right, see [20] and references therein) .

Our approach will allow one to study properties of trapped bosons, systems with a static or mobile impurities [21]. Also, it will be interesting to investigate three-dimensional bosonic bound clusters that appear in different branches of physics such as ${}^{4}\mathrm{He}$-clusters in condensed matter physics [22] and $\alpha$-clusters in nuclear physics [23, 24]. In these cases one might need to pick the basis carefully to reduce numerical effort. For instance, if the system is spherically symmetric then the basis should be chosen accordingly (cf. Ref. [14]).

With some modifications our method can be used to study other set-ups. In particular, we believe that it is possible to extend the method to bosonic systems without a condensate. To this end, one shall simply follow the steps presented above. First a reference state is used to normal order the operators. This reference state should describe an eigenstate of the Hamiltonian “well”, such that higher-body excitations are suppressed. As in the present work, the normal ordering provides one with means to truncate the differential equations, opening up the opportunity to approach $N$-body problems using a few-body machinery. Note that a suitable reference state in one-dimensional systems can be obtained by a linear superposition of weakly- and strongly-interacting states [25], providing one with a good starting point for this investigation.

We thank Achim Schwenk and Nikolaj Zinner for useful conversations. We thank Oleksandr Marchukov and Marcel Schmidt for their comments on the video abstract. A. G. V. gratefully acknowledges the support of the Humboldt Foundation. H.-W. H. was supported in part by the Deutsche Forschungsgemeinschaft through SFB 1245 and by the German Federal Ministry of Education and Research under contract 05P15RDFN1.

Appendix A Evaluation of commutators

To write down the flow equations, we need the commutators of the terms in ${\boldsymbol{\eta}}$ and ${\boldsymbol{H}}$. For the commutator of one-body operators and one- and two-body operators, we find:

[TABLE]

here we assume that $\beta_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}=\beta_{\alpha_{2}\alpha_{1}\alpha_{3}\alpha_{4}}=\beta_{\alpha_{1}\alpha_{2}\alpha_{4}\alpha_{3}}$, the same will be assumed for $\lambda_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}$. Note that with our definition of $\kappa=(N-1)/(2N)$ the element proportional to $\mathbf{I}$ in the second last row vanishes. For the commutator of the two-body operators, we find:

[TABLE]

The last term proportional to ${\boldsymbol{a}}_{i}^{\dagger}{\boldsymbol{a}}_{k}^{\dagger}{\boldsymbol{a}}_{b}^{\dagger}{\boldsymbol{a}}_{l}{\boldsymbol{a}}_{c}{\boldsymbol{a}}_{d}$ does not fit in our approximation scheme and should be truncated. Our implementation of this truncation is discussed in B.

Appendix B Truncation of the three body operator

To truncate the three-body operator, we assume that the number of particles in the lowest state is large, and thus the main contribution to the ground state energy is due to the piece of ${\boldsymbol{a}}_{i}^{\dagger}{\boldsymbol{a}}_{k}^{\dagger}{\boldsymbol{a}}_{b}^{\dagger}{\boldsymbol{a}}_{l}{\boldsymbol{a}}_{c}{\boldsymbol{a}}_{d}$ which contains at least one operator $a_{0}^{\dagger}$ and one operator $a_{0}$. Because of the presence of the condensate, these operators are then treated as numbers, i.e.,

[TABLE]

where

[TABLE]

Appendix C Convergence of the two body energy

The delta function interaction leads to a cusp in the wave function at zero separation of particles. This non-analyticity implies that accurate results for observables can be obtained only with a large number of plane wave states. We illustrate this statement by plotting the convergence of the ground state energy versus the number of the one-body basis states for the Lieb-Liniger model with just two particles, see figure 4. This plot shows that even in the two-body system the convergence with $n$ is very slow if $g$ is large. For large $n$ the convergence pattern in the figure can be well approximated by

[TABLE]

where $A$ is some constant that depends on $g$. Note that this convergence is faster than in a harmonic oscillator [26, 27] where it is described by $\sim 1/\sqrt{n}$. As is apparent from the discussion below this difference is connected to a slower growth of the energy with $n$ in a harmonic trap compared to a ring.

To understand this convergence pattern let us assume that we have diagonalized the matrix for some cutoff $n$, and obtained the energy $E_{n}$ and the wave function $\Psi_{n}$. Now let us see what happens when we diagonalize the Hamiltonian for $n+2$. The corresponding matrix includes the matrix for $n$ coupled to the rest via

[TABLE]

where at least one of the states $n_{1}$, $n_{2}$ was not included in the matrix for $n$. We have assumed that $n$ is large so $\Psi_{n}\simeq\Psi_{\infty}$. The function $\Psi_{\infty}(x_{1},x_{1})\equiv\bar{\Psi}$ is constant due to the rotational symmetry of the ring, and therefore we have

[TABLE]

Now using the second order correction from matrix perturbation theory we calculate the correction to $E_{n}$ due to the increase of the matrix size

[TABLE]

Summing the contributions for different $n$ up to infinity, this equation leads directly to (37). In general the leading order correction proportional to $n^{-\delta}$ is characteristic for delta function interactions and we can use it to obtain accurate results from the convergence pattern. We have observed that for larger number of particles the convergence behavior is also well described by (37).

Appendix D Effective Interaction

Another way to produce accurate results for the Lieb-Liniger model is to use some effective potential that reproduces low-energy properties of the system. For relevant studies of cold atomic systems see references [28, 29]. To introduce this potential, we first notice that all that we need to know about the interaction in our formalism is the following matrix element

[TABLE]

where $|n_{i}|\leq n_{max}$ and $n_{max}$ is the truncation parameter defined by $n$. Apparently such a matrix element also appears when we solve the Schrödinger equation in the ’relative’ coordinates

[TABLE]

by expanding the wave function $\psi$ in the plane wave basis, i.e., $\psi=\frac{1}{\sqrt{L}}\sum_{|n_{l}|\leq n_{max}}a_{l}e^{-2\pi in_{l}x/L}$ and solving the matrix equation

[TABLE]

Here $\mathbfcal{Q}$ is the matrix that contains eigenvectors as columns, ${\bf E}$ is the diagonal matrix that contains the eigenvalues, and $\mathbf{T}$ is the kinetic energy. Now we can turn the question around and find the potential that within our truncation space gives some specific matrices $\mathbf{E}$ and $\mathbfcal{Q}$. Such the potential then reads

[TABLE]

Let us now specify the desired low-energy properties. First of all, we fix the energies $E_{\alpha\alpha}$ to the $n$ lowest eigenenergies of the equation

[TABLE]

this choice means that in the two-body sector we always obtain correct energies. Next, we define the matrix $\mathbfcal{Q}$ as

[TABLE]

where the matrix ${\bf u}$ is an $n\times n$ matrix defined as $u_{\alpha_{1}\alpha_{2}}=\frac{1}{\sqrt{L}}\int_{0}^{L}e^{2\pi i\alpha_{1}x/L}\phi_{\alpha_{2}}(x)$. We see that if $n\to\infty$ then ${\bf u}^{T}{\bf u}\to 1$, and we have $\mathbfcal{Q}\to{\bf u}$. Therefore, the matrix ${\mathbfcal Q}$ is an orthogonal matrix that approximates the eigenstates and for $n\to\infty$ it reproduces the exact results.

The effective interaction shows faster convergence than the zero-range interaction, see figures 5 and 6, which depict our results for a few representative cases. By comparing the fitted values for the zero-range interaction and for the effective interaction we cross-check the two methods and insure accuracy of our results. The convergence pattern for the delta function potential is usually well described by 37. Note that we cannot directly apply the same line of arguments to find the convergence pattern for the effective interaction potential. Indeed, in this case the increase of the matrix size leads to a change of all matrix elements, and, therefore, standard perturbation theory cannot be used.

Appendix E Other generators

In section 2, we present examples of different generators that can be used to create the flow, see also [12, 13, 14, 15, 16]. In the main text, we illustrate our method using exclusively the operator (30) and leave other generators for future studies. Note that other ${\boldsymbol{\eta}}$ can be used directly in the derived flow equations (21)-(23) after the parameters of the generator (17) are specified. In this appendix, we briefly discuss the use of the White-type generator

[TABLE]

where

[TABLE]

This generator is similar to the one in (30) but it has additional energy denominators. As can be infered from section 2 for weak interactions this leads to the simultaneous decay of all couplings with $e^{-s}$.

Without truncation, the operators $\eta^{White}$ and $\eta$ in (30) define a unitary transformation and consequently lead to the exact energies. Our truncation scheme spoils this property, but it turns out that for the considered cases the results of the two generators are still very close to each other. We illustrate this statement in figure 7 for $N=4$ and $\gamma=\pi^{2}/2$. The correction $\delta e$ for this case accounts for about quarter of the SRG result meaning that the truncation procedure is no longer accurate, still the relative difference between the two results is a fraction of a percent. Therefore, for this problem these two generators can be used interchangeably.

Appendix F Dependence of $\delta e$ on $N$.

In our work we noticed that $\delta e$ from (31) increases with the number of particles $N$ for a fixed density $n$ and $\gamma$. This feature can be observed in figures 2 and 3 where for same values of $\gamma$ these corrections in the $N=4$ case are smaller than in the $N=15$ case. We report a similar behavior also in [21]. To understand this growth, let us first analyze the flow equations (21)-(23) with the generator (30) in the limit $gN\to 0$. In this case $S_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}\ll NS_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}\alpha_{5}\alpha_{6}}$, and, thus,

[TABLE]

We can estimate $\Gamma_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}$ using $f_{\alpha_{1}\alpha_{2}}(s)$ instead of $f_{\alpha_{1}\alpha_{2}}(0)$ in $S_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}$. If we do so, we find that $\Gamma$ is simply proportional to $g$, and, hence, $\delta e/e\sim\gamma^{4}N^{4}$. We see that in the limit $gN\to 0$ the correction grows very rapidly with $N$.

We are not able to provide a simple analytical analysis if the terms with $NS_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}\alpha_{5}\alpha_{6}}$ in (23) are large. Instead, we investigate this case numerically. To this end, we choose to work with $\gamma=0.1$. We find (see figure 8) that the ratio $\delta e/e$ increases with $N$, however, slower than $N^{4}$. Fitting suggests a much milder $\sim N^{2}$ scaling in this window of $N$.