Efficient construction of many-body Fock states having the lowest energies

TL;DR

This paper introduces a fast, simple algorithm for constructing many-body Fock bases from low-energy states, significantly improving the accuracy of exact diagonalization calculations for quantum systems.

Contribution

The authors present a novel algorithm that efficiently generates low-energy many-body Fock states, enhancing the accuracy of quantum many-body calculations.

Findings

The algorithm is insensitive to single-particle energy distributions.

It can be applied to bosonic and fermionic systems with arbitrary particle numbers.

Exact calculations using the new basis outperform standard approaches.

Abstract

To perform efficient many-body calculations in the framework of the exact diagonalization of the Hamiltonian one needs an appropriately tailored Fock basis built from the single-particle orbitals. The simplest way to compose the basis is to choose a finite set of single-particle wave functions and find all possible distributions of a given number of particles in these states. It is known, however, that this construction leads to very inaccurate results since it does not take into account different many-body states having the same energy on equal footing. Here we present a fast and surprisingly simple algorithm for generating the many-body Fock basis build from many-body Fock states having the lowest non-interacting energies. The algorithm is insensitive to details of the distribution of single-particle energies and it can be used for an arbitrary number of particles obeying bosonic or fermionic statistics. Moreover, it can be easily generalized to a larger number of components. Taking as a simple example the system of two ultra-cold bosons in an anharmonic trap, we show that exact calculations in the basis generated with the algorithm are substantially more accurate than calculations performed within the standard approach.

Equations

$$\hat{H}=\sum _{i=1}^N\left[-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial {\mathbf{r}}_i^2}+V({\mathbf{r}}_i)\right]+{\hat{H}}_{int}.$$

$$|\mathtt{i}⟩=|n_0{n}_1\dots ⟩,$$

$${\epsilon }_0$$

$${\epsilon }_k$$

$${\epsilon }_0$$

$${\epsilon }_k$$

$$\hat{H}=\sum _{i=1}^N\left[-\frac{1}{2}\frac{{\partial }^2}{\partial x_i^2}+\frac{1}{2}|x_i{|}^{2-\alpha }\right]+g\sum _{i>j}\delta (x_i-x_j),$$

Full text

1 Introduction

The exact diagonalization of the many-body Hamiltonian [1, 2, 3] is one of the simplest and straightforward methods of finding the ground state of the system of interacting quantum particles. It relies on a simple observation that having defined a finite set of ${\cal D}$ Fock states $|\mathtt{i}\rangle$ ($\mathtt{i}\in 1,\ldots,{\cal D}$) one can calculate all matrix elements of the Hamiltonian $\hat{H}$ of the system, $H_{ij}=\langle\mathtt{i}|\hat{H}|\mathtt{j}\rangle$, and numerically diagonalize corresponding matrix. As a result one obtains approximate decomposition of system’s eigenstates into defined Fock basis $|\Psi_{i}\rangle=\sum_{i}\alpha_{i}|\mathtt{i}\rangle$ and their eigenenergies $E_{i}$. Obviously, the accuracy of this straightforward method strongly dependents not only on the number of Fock states used for calculations but also on a particular method the states are chosen from the infinite set of all possible Fock states.

In typical physical scenario the total Hamiltonian of the many-body system can be divided to the sum of single-particle Hamiltonians describing dynamics of particles confined in an external potential $V(\boldsymbol{r})$ and the remaining part describing mutual interactions between particles, i.e., the many-body Hamiltonian of $N$ particles has a form

[TABLE]

Here we do not specify what is the form of the interaction Hamiltonian $\hat{H}_{int}$. In the most general case, it nontrivially depends on all particles’ positions $\{\boldsymbol{r}_{1},\ldots,\boldsymbol{r}_{N}\}$ and their derivatives. In the following, we assume that the single-particle part of the Hamiltonian is already diagonalized and all its eigenstates $\varphi_{k}(\boldsymbol{r})$ and their eigenenergies $\varepsilon_{k}$ (sorted along ascending order) are known. Having these states, it is very convenient to consider many-body Hilbert space of $N$ particles as spanned by Fock states build from these single-particle orbitals. It means that any Fock state can be written formally in the second quantization formalism as a sequence of numbers of particles occupying individual orbitals

[TABLE]

with $\sum_{k}n_{k}=N$. Of course, in the case of fermionic particles, an additional constraint has to be imposed assuring that for any $k$ occupation $n_{k}\in\{0,1\}$. In this way, quantum indistinguishability and statistics are taken into account. The non-interacting energy of this state is $E_{\mathtt{i}}=\sum_{k}n_{k}\varepsilon_{k}$. This notation immediately suggests one of the simplest methods of limiting the size of the Fock space and choosing only ${\cal D}$ of them as required by any numerical method based on diagonalization [4]. Namely, one limits the number of single-particle orbitals to some chosen cut-off ${\cal C}$, $k\in\{0,\ldots,{\cal C}\}$, and use them to build all possible Fock states for further calculations. In this approach, to increase the numerical accuracy of the many-body diagonalization, one slightly modifies a shape of single-particle orbitals [5] or increases the cut-off ${\cal C}$. Unfortunately, along with increasing cut-off ${\cal C}$ the number of Fock states grows exponentially (${\cal D}_{\mathtt{C}}=({\cal C}+N)!/({\cal C}!\,N!)$ and ${\cal D}_{\mathtt{C}}=({\cal C}+1)!/[({\cal C}-N+1)!\,N!]$ for bosons and fermions, respectively). In consequence, the method is naturally limited by numerical resources (mainly by available memory and computational time) and it is commonly viewed as very demanding from a computational perspective.

To limit exponential growth of the considered Hilbert space one can change the way the Fock basis is constructed. As discussed in the literature (see for example [6, 7, 8]), generating the Fock basis directly from the limited number of single-particle states is not the most efficient way of obtaining well-converged results, since then an energetic hierarchy in the Fock basis is completely neglected. Namely, this approach takes into account Fock states with relatively large energies fitting to small cut-offs $\cal C$ and neglects at the same time other Fock states with substantially lower energies but having larger single-particle excitations. As argued in [8], to increase the accuracy of the exact diagonalization without extending numerical efforts, one should select Fock states with the lowest non-interacting energy $E_{\mathtt{i}}=\sum_{k}n_{k}\varepsilon_{k}$ rather than states with the lowest excitations $k\leq{\cal C}$.

In general, selecting the Fock states $\{|\mathtt{i}\rangle\}$ having the lowest non-interacting energy is not a trivial task. Only in the case of equally distributed single-particle energies $\varepsilon_{k}-\varepsilon_{k-1}=\Delta$ (like in the one-dimensional parabolic confinement) it can be done quite easily since then the non-interacting energies can be represented uniquely by integer numbers. Consequently, the problem can be reduced to the task of finding different partitions of an integer (representing the non-interacting energy) into a fixed number of parts (representing individual particles) [9]. A detailed explanation is presented for example in [8]. Unfortunately, a simple generalization of this approach to cases when single-particle energies are not equally distant does not exist. Therefore, one needs to apply a direct selection of states with the lowest non-interacting energy from a basis of states having the lowest excitations (determined by a single-particle cut-off ${\cal C}$). This, however, is very inefficient and consuming huge amounts of computational facilities (memory and computation time). Consequently, its usefulness is strongly limited.

In this work, we present a simple numerical algorithm of generating the cropped Fock basis for $N$ particles (bosons as well as fermions) in the general case of any complete set of single-particle orbitals $\varphi_{k}(\boldsymbol{r})$ having ascending single-particle energies $\varepsilon_{0}\leq\varepsilon_{1}\leq\ldots$. As the input, the algorithm requires only one parameter, namely the maximal energy $\boldsymbol{E}$ of states in generated Fock basis. As a result one obtains a complete set of all ${\cal D}_{\mathtt{E}}$ Fock states having energy no larger than declared maximal energy $\boldsymbol{E}$. Importantly, an order of states in the generated set is well-known since automatically they are sorted lexicographically. Therefore, one can easily find not only $\mathtt{i}$-th Fock state but also perform reversed identification on demand and find the index $\mathtt{i}$ of the Fock state represented by its occupations. This property of the set is crucial for quick generation of a matrix representation of a many-body Hamiltonian matrix. A fundamental advantage of the algorithm presented is that Fock states are generated sequentially one-by-one in ${\cal O}({\cal D}_{\mathtt{E}})$ time, i.e., an improvement with respect to the case of a direct selection method is tremendous.

2 The algorithm

Counterintuitively, the best way to present the algorithm for generating the cropped Fock basis of states with the lowest non-interacting energy is to represent Fock states in the first quantization formalism, i.e., instead of the number of particles occupying single-particle orbitals one needs to remember single-particle states occupied by individual particles. It means that the Fock state $|\mathtt{i}\rangle$ is represented by a set of $N$ numbers $(i_{1},\ldots,i_{N})$ rather than a vector (2) encoding occupations. The relation between the two is straightforward. As an example, the non-interacting ground-state of $N=4$ bosons in both notations reads $|\mathtt{1}\rangle=|4\,0\,\ldots\rangle=(0,0,0,0)$ while in the case of $N=4$ fermions it has a form $|\mathtt{1}\rangle=|1\,1\,1\,1\,0\ldots\rangle=(0,1,2,3)$. Since particles are indistinguishable, it is understood that the states encoded in the first quantization are appropriately (anti-)symmetrized and only for simplicity they are represented by sets with ascending numbers, $i_{1}\leq i_{2}\leq\ldots\leq i_{N}$.

The idea of the algorithm is very similar for both quantum statistics. Therefore, let us first perform a short presentation for the case of $N$ bosons. Schematic flowchart of the algorithm (in the Fortran-like code) is presented in Fig. 1. The input requirements are: number of particles $N$, the single-particle energies $\varepsilon_{k}$ sorted along ascending order (stored in E(k)), and the maximal energy $\boldsymbol{E}$ (stored in EnergyMax). As an output one gets the series of Fock states $\{|\mathtt{i}\rangle\}$ represented by a set of numbers $(i_{1},\ldots,i_{N})$ (stored temporarily in State(j)) which are generated in lexicographical order. Example outputs of the algorithm for two different situations are presented in Table 1. Table 1a presents resulting Fock basis for the one-dimensional harmonic confinement having equally distant single-particle energy levels

[TABLE]

with the maximal energy $\boldsymbol{E}=6$. As it is seen, only twelve Fock states have this property (${\cal D}_{\mathtt{E}}=12$). At this point, one should note a huge reduction of the considered Hilbert space when compared to the standard approach of single-particle cut-off. It is clearly seen that in the case of the former method, to capture all states with non-interacting energy no larger than $\boldsymbol{E}=6$, one needs to use cut-off ${\cal C}=4$ (state $|\mathtt{5}\rangle$ has one particle excited to the single-particle state with energy $\varepsilon_{4}=4.5$). It means that the cropped Hilbert space generated with the standard cut-off approach would have a much larger dimension of ${\cal D}_{\mathtt{C}}=70$.

To show that the algorithm works also for other confinements in Table 1b we present the output for the same maximal energy $\boldsymbol{E}=6$, but for anharmonic confinement having decreasing energy distances between single-particle orbitals. As an example we set

[TABLE]

It is clearly seen that in this case, the gain from the algorithm is even larger. There are only sixteen Fock states (${\cal D}_{\mathtt{E}}=16$) having appropriate energy, while some of them have highly excited particles. For example, the state $|\mathtt{8}\rangle$ has one particle with energy $\varepsilon_{7}=92\,991/15\,625$. It means that in the standard cut-off method one needs to introduce ${\cal C}=7$ to capture all the states and the resulting size of the Hilbert space would be considerably larger ${\cal D}_{\mathtt{C}}=1\,980$.

In the case of fermionic statistics, any two particles cannot occupy the same single-particle orbital. Therefore, the algorithm needs to be modified to take into account this constraint. Fortunately, the changes are not substantial and they concern only the Initial Block and the Reset Block (marked in Fig.1 as INI and RES, respectively). In the Initial Block, the change is trivial and it is reduced to a simple change of the non-interacting ground state. In the Reset Block, particles having indexes larger than currently stored index should have consecutive occupations instead of occupations equal to the occupation of the indexed particle. For completeness, in Fig. 2 we present a schematic flowchart of the algorithm with these two changes directly incorporated.

Finally, let us emphasize that the algorithm presented can be very easily generalized to cases of multicomponent mixtures of particles with different statistics. One may think that the simplest way of such generalization is to generate Fock basis for each component separately and then build a multicomponent basis by a simple tensor composition. However, this approach is not correct since the resulting Fock basis does not contain states with the lowest non-interacting energy. It is especially crucial when the numbers of particles in different components are different or when different components have different single-particle energy spectra. To overcome this difficulty one should apply above algorithms to each component recursively taking care on the total non-interacting energy of generated Fock state, i.e., when the state of a given component is generated (instruction write(,) State in flowcharts) next algorithm for the following component should be executed. To assure that the non-interaction energy of the complete Fock state is bounded, the EnergyMax variable for a given component should be set as the maximal non-interacting energy $\boldsymbol{E}$ lowered by a sum of current energies of components generated already. The complete Fock state is accepted as a member of the Fock basis in the last component output.

3 Physical example

To show that the presented method of building the many-body Fock basis substantially increases efficiency of the exact diagonalization method let us consider the simplest model of $N=2$ ultra-cold bosons confined in a one-dimensional anharmonic trap and interacting via contact forces. In the dimensionless units the Hamiltonian of the system reads

[TABLE]

where $\alpha$ controls an anharmonicity of the trap (for $\alpha=0$ the standard harmonic confinement is restored) and $g$ is the interaction strength between particles. In the following example we set $\alpha=0.5$ and $g=0.1$.

In Fig. 3 we present the ground-state energy of the Hamiltonian (5) obtained with two complementary methods. First, with blue points, we display the ground-state energy obtained when the Fock basis is generated with the presented algorithm. Different points correspond to different maximal energies $\boldsymbol{E}$ and they are labeled by corresponding sizes of obtained Fock basis ${\cal D}_{\mathtt{E}}$. For comparison, using red horizontal lines, we present the ground-state energies obtained with the standard method of truncating the single-particle basis. Different lines correspond to different cut-offs ${\cal C}$ and they are additionally labeled by corresponding sizes of the cropped many-body Hilbert spaces ${\cal D}_{\mathtt{C}}$. It is clearly seen that along with increasing sizes of the Fock basis the ground-state energies decreases (ground-state energies obtained with any approximate method are always bounded from below by the exact ground-state energy of the system). Note that, the standard cut-off method is very inefficient when it is compared to diagonalization in the basis of energetically the lowest Fock states. For a given $\cal C$ the same ground-state energy may be obtained for much smaller Hilbert spaces provided that one selects Fock states systematically according to the non-interacting energy. For example, the ground-state energy obtained for ${\cal C}=26$ corresponding to the size of the Hamiltonian matrix ${\cal D}_{\mathtt{C}}=23\,751$ is almost the same as the energy obtained for almost four times smaller Hamiltonian matrix (${\cal D}_{\mathtt{E}}=6\,536$) calculated in the basis of Fock states with non-interacting energy not larger than $\boldsymbol{E}=26$.

4 Conclusions

Accuracy of the exact diagonalization of many-body Hamiltonians crucially depends on an appropriate selection of the Fock states. Typically, they are chosen from the set of eigenstates of the non-interacting system. Then, the efficiency of the method can be significantly improved if, instead of states having particles with the lowest excitations, the states with the lowest non-interacting energies are chosen. In this work, we present a very efficient and surprisingly simple algorithm generating all Fock states with the lowest non-interacting energy from a given set of single-particle orbitals. In this way, we give the base for further efficient calculations of a variety of many-body scenarios. In principle, the algorithm can be adapted to any many-body system for which a non-interacting Hamiltonian can be written as a sum of independent single-particle Hamiltonians. Therefore, it can be extremely useful in many different physical scenarios.

Acknowledgments

The authors would like to thank Przemysław Kościk for constructive comments. This work was supported by the (Polish) National Science Center through Grant No. 2016/22/E/ST2/00555.