Abstract

Highlights

1. Introduction

The solution of large, sparse eigenvalue problems is an important task in engineering, mathematics, and physics, particularly in the field of strongly correlated quantum many-body systems, such as the high-$T_C$ superconductors, and more recently cold atoms on optical lattices and coupled light-matter systems. The treatment of strong correlations between particles leads to exponentially large eigenvalue problems as a function of the system size.

To begin with, algebraic eigensolvers like the Lanczos and Arnoldi methods [1,2] play an important role in many-body physics as well as in many other research areas, mainly due to the fact that they are widely applicable, simple, effective and they yield numerically exact results. Implementations of the methods can be obtained from the Internet, e.g. the well-known ARPACK¹ routine which includes an implicitly restarting Arnoldi algorithm. On the downside, only comparably small systems can be treated by these techniques. For most problems in quantum many-body physics, the corresponding geometric sizes are much too small.

A couple of sophisticated numerical methods have been developed for such systems in recent years. The density matrix renormalization group (DMRG) [3] is a powerful algorithm for the determination of ground state properties of large systems, which are not algebraically feasible. Recently the approach has been embedded into a wider mathematical framework, the matrix product states (MPS, [4]). The approach is limited to 1D or quasi-1D systems. Another method is the cluster perturbation theory (CPT, [5]), which constructs approximations to the Greenʼs function of infinite systems in 1, 2, and 3 dimensions by exact treatment of small clusters and combining them by perturbation theory. An extension of CPT represents the variational cluster approach (VCA, [6]) which improves the results by introducing variational parameters. An approach without systematic errors is given by the family of quantum Monte-Carlo techniques (QMC, [7]). QMC simulations, which are based on high-dimensional random samples of some suitable probability-density function have a statistical error which declines with the sample size. They are applicable to fairly large systems and to finite temperatures. The drawback of these methods is the so-called sign-problem, that shows up especially in fermionic models and which makes certain models or parameter regimes inaccessible [8]. The methods discussed so far are particularly tailored for quantum many-body problems and cannot be applied easily, if at all, to matrices of other applications.

The methods for strongly correlated quantum many-body systems rely on the assumption that the model has only a few rather local degrees of freedom. Real ab-initio models for strongly correlated quantum many-body systems are out of reach in any case, and it is inevitable to reduce the complexity of the system upon describing the key physical properties by a few effective degrees of freedom, like in the multi-band Hubbard model [9]. In many cases, these models are still too complicated and cannot be solved reliably neither analytically nor numerically. For these cases it has been proven very useful to construct effective Hamiltonians like the Heisenberg- or tJ-model. They are obtained via the projection technique [9,10] upon integrating out such basis vectors which correspond to high energy excitations, or rather which have very large diagonal matrix elements in a suitable basis. Of course, it is desirable that the quantitative results of the effective model are close to those of the original model. But more than that, it is the qualitative generic physics of strongly correlated fermions at low energies, which one wants to understand. Since the original model itself is tailored for that purpose and is already a crude approximation of the underlying ab-initio Hamiltonian, it suffices to have an effective model that still includes the key physical ingredients in order to describe competing effects of strongly correlated quantum many-body systems. The standard projection technique defines the space of dynamical variables (basis vectors) which are of crucial importance for the low lying eigenvalues, which are marked by small interaction energy (small diagonal matrix elements). The residual dynamical variables (basis vectors) are treated in second order perturbation theory. If the model contains too complicated terms, such as the density assisted next-nearest neighbor hopping in the tJ-model, they are omitted. The resulting model has a significantly reduced configuration space and is solved by one of the above mentioned techniques.

Here we present a numerical scheme, that represents a threefold generalization of the projection technique: (a) it combines the two steps of the construction of the effective model and exact diagonalization, (b) no approximations are made as far as the effective model is concerned and (c) it allows a systemic inclusion of higher order terms up to the convergence to the exact result. The projection step is based on the Schur complement and results in a non-linear eigenvalue problem, which is solved exactly.

The numerical projection technique (NPT), that shall be discussed in this paper, has several advantages. First of all, it is not necessary to formulate an explicit effective Hamiltonian, which is not trivial in more complex models involving several bands and other degrees of freedom than just charge carriers [11]. And last but not least, NPT allows to systematically go to higher orders, which becomes necessary if the coupling strength is merely moderate.

NPT is applicable to other large-scale eigenvalue problems as well. The underlying matrix has to be sparse and the diagonal elements, after suitable spectral shift, have to fulfill the following criterion: a few diagonal elements are zero or small and the vast majority of the diagonal elements is (much) greater than the sum of the respective off-diagonal elements.

The method is demonstrated by application to two representative problems of strongly-correlated quantum many-body physics, the spinless fermion model with nearest neighbor repulsion and the drosophila of solid state physicists: the Hubbard model with a local Coulomb repulsion. These models are used for benchmark purposes only and it is not the goal of the present paper to provide a detailed discussion of the fascinating physics described by these models.

The outline of the paper is as follows. In the second section we present the two benchmark models. The numerical projection technique is introduced in detail in the third section and analyzed in the following section. A further extension of the algorithm to improve the accuracy will be presented in Section 5, Section 6 will show results for the Hubbard model and concluding remarks will be given in the last section.

2. Model systems and basic idea

The approach that we present below is generally applicable if the matrix, for which the lowest eigenvalues shall be determined, can be split into parts of increasing diagonal dominance. What this means in detail will be clarified using two examples of the realm of many-body physics, the spinless fermion model with strong nearest neighbor interaction and the Hubbard model for fermions. Both models are tailored to study effects of strong correlations of electronic systems, such as the Mott-insulators [12], high temperature superconductors [13], manganites [14], just to name a few of the very many novel materials with fascinating many-body effects. The Hamiltonian of the spinless fermion model reads

where ${\hat{a}}_i^{†}$ (${\hat{a}}_i$) denotes the creation (annihilation) operator for fermions at site i. These operators have the common fermionic anti-commutator relations (for details see [15]). The operator ${\hat{n}}_i={\hat{a}}_i^{†}{\hat{a}}_i$ is the particle number operator for site i. The bracket $\langle ij\rangle$ indicates that the sum is restricted to nearest-neighbor sites $x_i$ and $x_j$. There is a total of L lattice sites $x_i$, which are placed on a one-dimensional regular lattice in the present work.

In the case of the Hubbard model, the spin degree of freedom is also taken into account and the Hamiltonian reads

The creation (annihilation) operators ${\hat{a}}_{i\sigma }^{†}$ (${\hat{a}}_{i\sigma }$) obtain a spin index σ with two possible orientations ↑, ↓. The operator ${\hat{n}}_{i\sigma }={\hat{a}}_{i\sigma }^{†}{\hat{a}}_{i\sigma }$ is the particle number operator for spin σ at site i and the operator for the particle density at site i is given by ${\hat{n}}_i={\sum }_{\sigma }{\hat{n}}_{i\sigma }$.

The physical behavior of both models is determined by the relative strength of the hopping parameter t and the two interaction parameters U, which stands for the on-site Coulomb repulsion between fermions of opposite spin, and V, which represents the repulsive nearest-neighbor interaction. In order to simplify the discussion we denote both parameters by V. In both cases the most interesting case is that of (almost) half-filling. That is, the number of particles in the system is half the maximum capacity, which is L in the spinless fermion model and 2L in the Hubbard model. Both models are particularly interesting in the strong coupling case, i.e. for $V>|t|$.

For both models we use the occupation number basis in real space in which the interaction part is diagonal and has a value $N_{\kappa }V$, where $N_{\kappa }$ is either the number of occupied nearest-neighbor sites ($N_{nn}$) in case of the spinless fermion model, or the number of double-occupancies ($N_d$) in case of the Hubbard model. In strong coupling, states with increasing values of $N_{\kappa }$ are decreasingly important. In the projection technique [9,10] effective strong coupling models are derived, in which the configuration space is restricted to the sector $N_{\kappa }=0$, i.e. $N_{nn}=0$ ($N_d=0$) and the influence of the higher sectors are taken into account up to second order in $t/V$. A prominent example is the spin-1/2 antiferromagnetic Heisenberg model which is obtained in the half-filled case of the Hubbard model or the tJ-model, which is the generalization away from half-filling. Already in the tJ-model terms are neglected which belong to the same order in $t/V$ and to the same $N_d$-sector. In models with more bands and degrees of freedom, the derivation of an effective strong coupling model can be very demanding, like e.g. in the spin-orbital model for the manganites [11].

In this paper, we exploit this idea numerically and generalize it in such a way, that the influence of high sectors is taken into account recursively without the need of leaving terms out that in standard projection technique would complicate the resulting effective Hamiltonian. For obvious reasons we will refer to this approach as numerical projection technique (NPT).

3. Numerical projection technique

In order to exploit the NPT, we reorder the basis vectors according to $N_{\kappa }$. The corresponding Hamiltonian matrix has a natural block structure (see Fig. 1) corresponding to the sectors with $N_{\kappa }=0,1,2,\dots .$

Open in a separate window

Fig. 1

Illustration of the Hamilton matrix of a V-ordered occupation number basis. The blue parts contain the hopping terms inside a partition and the potential on the diagonal and the green parts contain the hopping terms between the partitions. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Note, that the hopping of a particle can only change the number $N_{\kappa }$ by one. Therefore, only neighboring sectors are coupled and the matrix has a tridiagonal block structure.

3.4. Truncation of the first partition

The advantage of the numerical projection technique so far is the fact that (non-linear) eigenvalue problems have to be solved for matrices of a size given by the first partition $p_1$, which is obviously smaller then the original size (see Fig. 2).

Open in a separate window

Fig. 2

Truncation procedure for the scheme presented in the text. Each partition is consecutively projected onto the original effective matrix and diagonalized.

It may happen that the first partition, even if it covers merely the first sector, is still very large. This is the case in the Hubbard model, where already the lowest sector for $N_{\kappa }=0$ yields $p_1=\left(\begin{array}{c}\hfill L\hfill \\ \hfill N_{\uparrow }\hfill \end{array}\right)$ $\left(\begin{array}{c}\hfill L-N_{\uparrow }\hfill \\ \hfill N_{\downarrow }\hfill \end{array}\right)$, which at half-filling becomes $p_1=\left(\begin{array}{c}\hfill L\hfill \\ \hfill L/2\hfill \end{array}\right)$. Next we will try and reduce the size of the first partition. Firstly, we can exploit translational invariance and solve the eigenvalue problem in the individual k-space sectors. We can furthermore restrict the number of vectors, which are retained in the first partition, based on some suitable criterion as to their importance for the lowest eigenvalue. We will use only the most obvious and simple criterion, namely to keep those vectors, which belong to the lowest eigenvalues $D_2$ of $A_1$. More sophisticated criteria as in DMRG are also conceivable.

We will use a predefined value ${\tilde{p}}_1⩽p_1$ for the truncation size of the first partition. That means, we solve the eigenvalue problem of $A_1$ of size $p_1\times p_1$ and retain only those ${\tilde{p}}_1$ eigenvectors in $V_1$ which correspond to the lowest eigenvalues. According to the recursion procedure, in the following recursion steps the diagonal matrices entering the left-upper block $D_n$ will all have the size ${\tilde{p}}_1$ and henceforth only (non-linear) eigenvalue problems of this size occur.

4. Numerical experiments

We begin with the spinless fermion model with periodic boundary condition (pbc) at half-filling. Due to the conservation of the total momentum k only those basis vectors with the same k need to be taken into account. For the spinless fermion model at half-filling the ground state is degenerate and has $k=\pm \frac{\pi }{2}$. We will therefore only take basis vectors into account, that have a total momentum $k=\pm \frac{\pi }{2}$.² The selection of the total momentum leads to a significant reduction of the computational effort. On the one hand, the size $p_1$ of the initial partition is reduced roughly by L and on the other hand it can be avoided that in the truncation step vectors are retained, that are irrelevant for the ground state of the original system.

To begin with, we choose as first partition the first sector with $N_{\kappa }=0$, which consists of merely 2 states with particles on every other site (charge density wave), i.e. ${\tilde{p}}_1=p_1=2$. These two vectors are the real valued linear combinations of the $\pm \frac{\pi }{2}$ basis vectors.

Table 1 exhibits a comparison of the run time for various system sizes. It can be seen that NPT is significantly faster than the highly optimized ARPACK routine (Arnoldi implementation) for increasing system size. We observe that even for ${10}^6\times {10}^6$-matrices the CPU time is less than 1 sec and more than a hundred times faster than the highly optimized ARPACK code.

4.1. Dependence on the size of the initial partition

We observe that the convergence of the eigenvalue levels off at a moderate relative accuracy of 1 percent when the initial partitions consist only of the first sector ($N_{\kappa }=0$, see Fig. 3). The reason is due to the fact, that the vectors that are omitted during the recursion procedure, are relevant for a higher accuracy.

Open in a separate window

Fig. 3

Convergence of the lowest eigenvalue with the largest sector ($N_{\kappa }^{\max }$) included in the recursion, depending on the initial partition (e.g. green squares: the initial partition consists only of the first sector $N_{\kappa }=0$). In all cases, ${\tilde{p}}_1=2$ was used, i.e. the 2 lowest eigenvectors of the initial partition where kept. The blue solid line indicates the exact values corresponding to the various subspaces spanned by the vectors of the sectors $N_{\kappa }=0,\dots ,N_{\kappa }^{\max }$, respectively. The results are compared with the exact eigenvalue $E_0$ of the original matrix on a logarithmic scale. The system parameters are L = 20,N = 10,|V/t| = 10. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

One way to increase the accuracy is to choose a greater initial partition upon including all sectors up to $N_{\kappa }^{\mathrm{ini}}$, but still using a truncation size of the first partition ${\tilde{p}}_1=2$. See the curve with turquoise upward triangles in Fig. 3 for results with initial partition containing sectors up to $N_{\kappa }=2$. The relative accuracy increases by one order of magnitude but is still restricted. In all curves we observe that the result is converged when two more sectors are included beyond the ones present in the initial partition, which is obviously a consequence of the small truncation size ${\tilde{p}}_1=2$.

5. Three-partition projection

As discussed in the introduction, the central goal of the projection techniques is to describe the low-lying generic physics of the original quantum many-body system. We have seen, that NPT as opposed to standard projection techniques resulting in effective models, even yields quantitatively good agreement with the result of the original model. If still higher accuracy is the ultimate goal, then one can do even better then what has been discussed so far. It is possible to include in each recursion step two sectors at a time with little more CPU effort. To this end we consider a matrix which has $3\times 3$ tridiagonal block structure, that results, if we add in each recursion (projection) step two additional sectors to the first partition. The matrix has the structure

Here, an equivalent relation to Eq. (3) can be obtained:

The contribution of the matrix C shall be expanded up to the linear inverse term:

Again, the expression is expanded in terms of the kinetic part:

Using this extension, the numerical result can be improved by orders of magnitude as is shown in Fig. 4. The results are for the spinless fermion model at strong coupling $|V/t|=10$ and half-filling. In both case the initial partition is given by the lowest sector without further truncation, i.e. $p_1={\tilde{p}}_1=2$.

Open in a separate window

Fig. 4

Comparison of the two-partition and three-partition approach for two different system sizes. The convergence of the relative uncertainty of the lowest eigenvalue is shown versus the maximum number of sectors included.

6. Results for the Hubbard model

Next we present some results for the Hubbard model. Standard projection technique leads in the case of half-filling to the Heisenberg model and away from half-filling to the tJ-model [9]. These models correspond in NPT to the situation that the initial partition is the first sector with $N_{\kappa }=0$ (no double occupancies) and only the second sector is mixed in the projection step. On top of that, the non-linear eigenvalue problem is replaced by the linear one with $\lambda =0$ in Eq. (2). In Table 3 the lowest eigenvalues for different system sizes and particle numbers are given. In the half-filled case the NPT results are compared with the exact results for the corresponding Heisenberg model and away from half-filling with those of the tJ-model. In NPT we use the first sector ($N_{\kappa }=0$) as initial partition, restricted to $k=0$. For the truncation size in the projection steps we used ${\tilde{p}}_1=2$. The NPT recursion steps are stopped at an absolute accuracy ${\epsilon }_{\lambda }={10}^{-4}$.

We compare NPT results and those of the effective models with the exact lowest eigenvalue of the Hubbard model. We see that even in the strong coupling case ($|U/t|=10$), investigated here, NPT yields quantitatively better results than the effective models.

7. Conclusions

The numerical projection technique (NPT), presented in this paper, is a numerical generalization of standard projection technique, routinely used to derive effective models for qualitative description of the generic low-energy physics of strongly correlated systems. The standard projection technique corresponds to an approximate evaluation of second order perturbation theory in the inverse interaction parameter.

NPT is on the one hand a systematic application of the projection techniques to any many-body Hamiltonian, irrespective of the complexity of the original model, also in cases where the analytic approach would lead to unmanageable effective models. On the other hand it is also applicable for moderate coupling parameters, where second order perturbation theory is not reliable enough.

We have demonstrated that NPT yields fairly accurate results with little CPU-effort as compared to state-of-the-art exact diagonalization. Furthermore, by solving the non-linear eigenvalue problem and including higher order terms the approach yields significantly better results than traditional projection technique.

Footnotes

References