# Extensive v2DM study of the one-dimensional Hubbard model for large lattice sizes: Exploiting translational invariance and parity

###### Abstract

Using variational density matrix optimization with two- and three-index conditions we study the one-dimensional Hubbard model with periodic boundary conditions at various filling factors. Special attention is directed to the full exploitation of the available symmetries, more specifically the combination of translational invariance and space-inversion parity, which allows for the study of large lattice sizes. We compare the computational scaling of three different semidefinite programming algorithms with increasing lattice size, and find the boundary point method to be the most suited for this type of problem. Several physical properties, such as the two-particle correlation functions, are extracted to check the physical content of the variationally determined density matrix. It is found that the three-index conditions are needed to correctly describe the full phase diagram of the Hubbard model. We also show that even in the case of half filling, where the ground-state energy is close to the exact value, other properties such as the spin-correlation function can be flawed.

## I Introduction

The reduced density matrix makes its first appearance in the work of Dirac, in which the single-particle density matrix (1DM) is used in the description of Hartree-Fock theory Dirac (1930). Husimi Husimi (1940) was the first to note that, for a system of identical particles interacting only in a pairwise manner, the energy can be expressed exactly as a function of the 2DM. This becomes very clear in second-quantized notation (see e.g. Fetter and Walecka (1971); Dickhoff and Van Neck (2008)), where a system of identical particles interacting pairwise is described by a general Hamiltonian:

(1) |

The expectation value for the energy of any ensemble of -particle wave functions with positive weights , can then be expressed as a function of the 2DM alone:

(2) |

in which we have introduced the 2DM:

(3) |

and the reduced two-particle Hamiltonian,

(4) |

The idea to use the 2DM as a variable in a variational scheme was first published in literature by Löwdin in his groundbreaking article Löwdin (1955), but even earlier, in 1951, John Coleman tried a practical variational calculation on Lithium. To his surprise, the energy he obtained was far too low, after which he realized the variation was performed over too large a class of 2DM’s Coleman (2000). Independently and unaware of the work by Löwdin and Coleman, Joseph Mayer Mayer (1955) used the 2DM in a study of the electron gas. In a reply to Mayer’s paper, Tredgold Tredgold (1957) pointed out the unphysical nature of the results, and suggested that additional conditions on the density matrix are needed to improve on them.

These results led Coleman, in his seminal review paper Coleman (1963), to formulate the -representability problem. This is the problem of finding the necessary and sufficient conditions which a reduced density matrix has to fulfil to be derivable from a statistical ensemble of physical wave functions, i.e. expressible as in Eq. (3). In this paper he also derived the necessary and sufficient conditions for ensemble -representability of the 1DM, and some bounds on the eigenvalues of the 2DM. A big step forward was the derivation of the and matrix non-negativity conditions by Garrod and Percus Garrod and Percus (1964). These were practical constraints, which allowed for a computational treatment of the problem. The first numerical calculation using these conditions on the Beryllium atom Garrod and Fusco (1976); Garrod et al. (1975) was very encouraging, as the results were highly accurate. It turned out, however, that Beryllium, due to its simple electronic stucture, is a special case where these conditions perform very well. A subsequent study showed that these conditions do not work well at all for nuclei Mihailovic and Rosina (1974); Rosina and Garrod (1975). This disappointing result, together with the computational complexity of the problem, caused activity in the field to diminish for the next 25 years. The change came with the development of a new numerical technique, called semidefinite programming, which turned out to be very suited for the determination of the 2DM under matrix non-negativity constraints. Maho Nakata et al. Nakata et al. (2001) were the first to use a standard semidefinite programming package to calculate the ground-state energies of some atoms and molecules, and obtained quite accurate results. He was quickly followed by the extensive work of Mazziotti Mazziotti (2002). These results reinvigorated interest in the method, and sparked of a lot of developments. New -representability conditions were introduced, e.g. the three-index conditions, set forth by Zhao et al. Zhao et al. (2004), which led to mHartree accuracy Hammond and Mazziotti (2005); Nakata et al. (2008); Mazziotti (2005); Gidofalvi and Mazziotti (2006); Mazziotti (2007); Braams et al. (2007) for molecules near equilibrium geometries.

In recent years interest in the method has been growing, as the variational determination of the 2DM results in a lower bound, which is highly complementary to the upper bound obtained in variational approaches based on a wave-function ansatz. In addition, the method is essentially non-perturbative in nature, and has a completely different structure quite unrelated to other many-body techniques. A lot of activity has been devoted to the search for new -representability conditions, which improve the result in a computationally cheap way Van Neck and Ayers (2007); Verstichel et al. (2010). There have also been efforts to improve the semidefinite programming algorithms by adapting them to the specific problem of density matrix optimization Mazziotti (2004, 2011); Verstichel et al. (2011), allowing the study larger systems.

The one-dimensional Hubbard model Hubbard (1963) is the simplest model possessing non-trivial correlations present in a solid. The Hamiltonian reads:

(5) |

It pertains to a one-dimensional lattice, with sites labeled by . Periodic boundary conditions (PBC) are assumed.

The complexity of this seemingly simple schematic Hamiltonian lies in the competition between the first term, called the hopping term, which delocalizes the electrons, and the second on-site repulsion term which is diagonal in the site basis. In Figure 1 a graphic representation of the two terms in the Hubbard Hamiltonian is shown.

In this article the one-dimensional Hubbard is studied using the standard two-index and three-index non-negativity conditions. In the next Section we discuss how the huge amount of symmetry which is present in the model can be exploited to get a significant speedup of the programs. More specifically, it is shown how translational invariance and space-inversion parity can be combined. In the subsequent Section this symmetry is used to compare the computational scaling of different semidefinite programming algorithms with increasing lattice sites. In the final Section we show the v2DM results using two- and three-index constraints for the one-dimensional Hubbard model on a 20- and 50-site lattice with different filling factors. We have not only computed the ground-state energy, but also compared the spin and charge two-particle correlations functions with Quantum Monte Carlo Sorella et al. (1990) and Bethe ansatz Ogata and Shiba (1990) results, to check the validity of the variationally obtained 2DM.

## Ii Exploiting symmetry in the Hubbard model

As there is no preferred direction in spin space in Eq. (5), spin symmetry can be exploited and all the results from Verstichel et al. (2009) are taken over. However, far more symmetries are present in the model, which allow for a huge reduction of the matrix dimensions involved.

### ii.1 Translational invariance

When periodic boundary conditions are assumed (as in Figure 1) the Hamiltonian is invariant under translations along the lattice. This is an abelian symmetry for which it is easy to block diagonalize the 2DM and its matrix maps, as the correct basis transformation in single-particle space automatically block diagonalizes all matrices on higher-order particle space. Translational invariance is exploited by Fourier transforming the site basis to quasi-momentum space:

(6) |

where is the lattice size, and takes on the values for . The kinetic or hopping part of the Hamiltonian becomes diagonal in this basis:

(7) |

from which it follows that in the non-interacting () ground state, the electrons occupy the states with momenta lower than the Fermi level.

The eigenstates of the Hubbard Hamiltonian have a good total quasi-momentum . The 2DM for these states, expressed in the quasi-momentum single-particle basis:

(8) |

where

(9) |

is automatically block diagonal, because the only non-zero matrix elements in Eq. (8) are those which conserve momentum: ^{1}^{1}1Here % signifies the modulo operator. This means we have blocks for every , with two-particle states that satisfy , and a block dimension that scales linearly with lattice size .

The spin-symmetric matrix constraints simplify considerably by including translational invariance, because the 1DM is automatically diagonal in the quasi-momentum basis:

(10) |

The translationally invariant 1DM can be derived from the 2DM as:

(11) |

As an example the translationally invariant form of the condition is shown. There is a slight complication because the correct annihilation or hole operator is given by:

(12) |

Using this operator the translationally invariant map becomes:

(13) |

where and with the particle-hole operator defined by:

(14) |

The map can be expressed as a function of the 2DM:

(15) |

from which one can see that the blocks in the matrix with , correspond to the blocks with in the 2DM, as expected for a particle-hole transformed quantity.

### ii.2 Parity

The Hubbard model with periodic boundary conditions (PBC) has more symmetries than spin and translational invariance, one of them being parity. Parity follows from the symmetry under the inversion of space, i.e. . One can readily appreciate from Figure 1 that this symmetry is present here, i.e. the Hamiltonian is invariant for the inversion of the site index . From the Fourier transform in Eq. (6) it can be seen that the effect of this operator on a momentum state is to transform into . However, this operation does not commute with translation, which means that the eigenstates of the Hamiltonian cannot have good momentum and parity at the same time. From now on we assume that the number of lattice sites is even.

As a consequence, if the ground state has momentum or then it is doubly degenerate, forming a doublet with . In what follows we will use this degeneracy to exploit both translational invariance and parity to reduce the dimensions of the matrices involved in the program. The following considerations are valid for every . We start with the simplest case, the 1DM.

Single-particle space is built out of different momentum states having up or down spin, , for and . If we transform to a basis with good parity, momentum is no longer a good quantum number:

(16) |

Two states, and ^{2}^{2}2Note that we use for both the parity quantum number, as for the transcendent number; in what follows it is always clear from the context what interpretation should be given to .
are mapped on themselves, and only positive parity states can be formed with these momenta. They have norms . For the other states, , both positive and negative parity combinations can be constructed, with norm . To take advantage of both symmetries at the same time, we define the 1DM using an ensemble of the and states:

(17) | ||||

(18) |

in which is a regular translationally invariant 1DM as defined in Eq. (10). Because of parity symmetry we have,

(19) |

from which it follows that Eq. (18) is diagonal in . If we now define:

(20) |

the translationally invariant 1DM with good parity reduces to:

(21) |

which can be seen to be independent of parity. The 1DM is still diagonal in , but less elements have to be stored, since for there is a degeneracy in parity.

In contrast to parity in atomic systems (see e.g. Verstichel et al. (2009)), the parity of a two-particle state for translationally invariant systems is not the product of the parities of the single-particle states building up the two-particle state. Instead, the parity is inherently a two-particle property, and the single-particle states building up the two-particle states have no good parity:

(22) |

with

(23) |

In general the 2DM is defined using a ensemble:

(24) |

where

(25) |

and with defined as in Eq. (9). Because of this -ensemble definition and the fact that parity symmetry implies that,

(26) |

one sees that the 2DM becomes diagonal in two-particle parity. As was the case for the 1DM, the and are mapped on themselves, but because the single-particle momenta change, both positive and negative parity combinations can now be formed. Let us take a look at the different possibilities:

#### :

for the single-particle indices in Eq. (22) have to satisfy:

(27) |

This means that states can be written as:

(28) |

in which the second term is equal to the first but with exchanged single-particle indices. From previous discussions we know that the symmetry of the two-particle state under the exchange of the single-particle indices depends on the two-particle spin, i.e.

(29) |

One can see from Eq. (28) that for only the positive parity states, and for only negative parity states remain. The norm is given . In the case, the definition of the parity-symmetric 2DM as a function of the regular translationally invariant 2DM then reduces to:

(30) |

#### :

for the first and second term in Eq. (22) consist of different single-particle indices , implying that both positive and negative parity combinations can be constructed, with norm . As shown for the 1DM, the ensemble makes the 2DM diagonal in, and independent of, parity. Hence every block is twofold degenerate. Since there are only two terms remaining in the definition of the parity-symmetric 2DM:

(31) |

#### :

Finally, for this block again equals . In this case there is always one state that is mapped on itself, and for which only a positive parity combination can be formed, i.e. , with norm . For all the other states in this block both positive and negative parity combinations can be formed, with norm . Because of the ensemble, the 2DM falls apart in a positive and negative parity block, and since , four terms remain in the definition of the 2DM:

(32) |

We observe from Eq. (32) that the original block reduces to a positive and negative parity block, for both the and part. Also note that there is no degeneracy between the positive and negative parity block!

Similar considerations hold for the matrix constraints, as an example the explicit case of the condition is given.

The parity-symmetric form of a particle-hole state is defined as:

(33) |

in which the hole operator is defined as in Eq. (12). Using this parity-symmetric particle-hole operator the map is defined in a ensemble, which once again renders the matrix diagonal in particle-hole parity. The particle-hole states can be divided into two classes, on the one hand , and on the other hand those states which are mapped on a different momentum. For simplicity we first consider this last class.

#### :

in this case one can construct both positive and negative parity combinations, with norm . The resulting matrix contains only two terms, because of momentum conservion, and is independent of particle-hole parity, so every block is twofold degenerate:

(34) |

This implies the following expression of the map as a function of the parity-symmetric 2DM:

(35) |

#### and :

both the and blocks are mapped on themselves. For the action of the parity operator is again to exchange the single-particle momenta, but in contrast with the two-particle case, there is no symmetry between the particle and the hole index. As a consequence positive and negative parity combinations for both and can be constructed, with norms . There are a few exceptions however: for , the states with and , and for the states with , and , are mapped on themselves and only occur in the positive parity block, with norm . The general form of the parity-symmetric map when , as a function of the regular translationally invariant is:

(36) |

In this case the expression of as a function of the 2DM is a bit more complicated:

(37) |

## Iii Computational performance

In v2DM we want to optimize the energy by varying a matrix, the 2DM, under the constraints that it has the right particle number, and that some linear matrix maps of the 2DM are positive semidefinite, i.e.

(38) | |||||

u.c.t. | (39) |

Here the are a collection of matrix non-negativity constraints to be applied. It is well known that this problem can be reformulated as a semidefinite program. There is a vast literature on this subject and many different algorithms exist to solve this type of problem. Because of the large amount of symmetry present in the Hubbard model, there is a huge block diagonalization of the 2DM and the matrix constraints. This leads to a significant reduction of the computational cost of a basic matrix computation. This feature has allowed us to go up to large lattice sizes and compare the computational scaling of different algorithms. We have implemented three different semidefinite programming algorithms. Two are so-called interior point methods (a dual-only potential reduction algorithm, Verstichel et al. (2009) and a primal-dual interior point algorithm Verstichel et al. (2011)), where the 2DM is optimized from within the -representable region. In the third one (a boundary point method Mazziotti (2011)) the 2DM is not required to be -representable during the optimization. In this Section we compare the computational scaling of these methods for the one-dimensional Hubbard model with and using conditions. For details about the implementation of the different algorithms, see the cited references.

All the methods have the same basic computational scaling behaviour, being (with the dimension of single-particle Hilbert space) required for multiplying, inverting or diagonalizing a matrix. In Figure 2 the number of these operations needed to converge to the optimum is plotted as a function of lattice size. The interior point methods both have to solve a linear system of size , so it is not surprising that the scaling, on top of the basic matrix computations, of these methods is . More surprising is that there seems to be no, or a very limited, scaling for the boundary point method. The number of iterations required remains around 3000, irrespective of the size of the system. It must be stressed that this is a result limited to the one-dimensional hubbard model, and cannot be extrapolated to other systems, like molecules, where the convergence properties of the boundary point method can be completely different. One reason for the succes of the boundary point method applied to the Hubbard model is the amount of symmetry present in the system. The boundary point method is designed for problems with a huge amount of dual variables or primal constraints. For most physical systems the dimensions of the matrices involved are already unfeasibly large before the boundary point method would becomes advantageous. The Hubbard model, however, contains many symmetries, implying that the matrix dimensions are considerably reduced, and the number of dual variables can get very large before the matrix computations involved become unfeasible. In this case, the domain where the boundary point method is advantageous is actually reached.

## Iv Results

In this Section we present and discuss the results of v2DM calculations, taking advantage of all the symmetries, on a 50-site lattice with the conditions, and on a 20-site lattice with the () conditions. The Hubbard model has been studied before using the v2DM method, see e.g. Hammond and Mazziotti (2006); Shenvi and Izmaylov (2010); Nakata et al. (2008); Verstichel et al. (2012), but only the energy was considered, and this for relatively small lattice sizes (up to ). In this paper we study different filling factors, and extract various properties like the ground-state energy and two-particle correlation functions in order to assess the quality of the variationally obtained 2DM. The v2DM results discussed in this Section were all obtained using the primal-dual predictor corrector semidefinite programming algorithm Verstichel et al. (2011). Although the one-dimensional Hubbard model can be solved exactly using the Bethe ansatz Bethe (1931); Lieb and Wu (1968); Essler et al. (1991, 2005), it is hard to extract information about the solution for finite systems. For the calculations on a 20-site lattice, we compare the data with the quasi-exact results obtained through a variational Matrix Product State (MPS) algorithm Schollwöck (2011); Verstraete et al. (2008); Chan and Zgid (2009), written by co-worker Sebastian Wouters Wouters et al. (2012). For the 50-site lattice, however, this is no longer computationally feasible. At half filling a simplification in the Bethe-ansatz equations occurs, which allows to calculate the ground-state energy of finite systems by solving a set of non-linear equations (Lieb-Wu)Lieb and Wu (2003). At other fillings no data is available for comparison.

### iv.1 Ground-state energy

In Fig. 3 the ground-state energy per particle of the one-dimensional Hubbard model is plotted as a function of the on-site repulsion (the hopping parameter will always be taken equal to unity). In the top figure the v2DM results for the 20-site lattice are shown for three different fillings, 12 particles (), 16 particles () and half filling. These were calculated using both the and the conditions, and are compared to the quasi-exact variational MPS results. In the bottom figure the v2DM results for the 50-site lattice are shown for the same fillings (i.e. 30 particles (), 40 particles () and half filling). For the 50-site lattice it was only possible to perform the calculations using the conditions, and compare to the exact solution obtained by solving the Lieb-Wu equations for the half-filled lattice Lieb and Wu (2003).

One interesting thing to notice is that the energy per particle for the 20-site lattice and the 50-site lattice, at the same filling, are very similar. This is due to the periodic boundary conditions which make the results converge quite rapidly for increasing lattice size , implying that one can already extract relevant results for the thermodynamic limit by studying relatively small lattices. This fast convergence can be clearly seen in Fig. 4, where we plotted the energy per particle of a Hubbard model with at half filling, as a function of the lattice size . This result seems to indicate that the method is more or less size extensive for the Hubbard model, which is surprising since in general v2DM is not size extensize Nakata and Yasuda (2009); van Aggelen et al. (2009).

Another thing to remark in Fig. 3 is that for the 20-site lattice, the difference between and is rather small for the half-filled lattice, but larger for the other fillings, and that the difference gets larger when increases. For the 50-site lattice we see that the result agrees nicely with the solution of the Lieb-Wu equations. For the other fillings no reference data are available. There are, however, two limits of the model that are exactly solvable. The first limit is the rather trivial case of no interaction, i.e. , for which the solution has already been given in Eq. (7). The Hamiltonian reduces to a single-particle operator, which means this limit is already described correctly by including the and conditions alone. The other exactly solvable limit is when . In this limit the physics of the model decouples into two independent parts, one describing the spin of the system, and the other the movement of the particles (this is called spin-charge separation Ogata and Shiba (1990)). This decoupling shows up in the Bethe-ansatz wave function: the charge degrees of freedom are described by a Slater determinant of spinless fermions, whereas the spin degrees of freedom become equivalent to a spin- Heisenberg model. The single-particle energy spectrum changes slightly compared to Eq. (7) because the boundary conditions for spinless fermions are periodic/antiperiodic if is even/odd Ogata and Shiba (1990); Krivnov et al. (1990):

(40) |

When the lattice is half-filled all the single-particle states are occupied, and the total energy sums up to zero, which is correctly described by the results in Fig. 3. Away from half filling, however, the energy has a finite limit which can be calculated using Eq. (40). From the figure we can see that the conditions do not suffice to correctly describe the large- limit. Only when the conditions are added, the results converge to the right limit. Calculations at very large values of have been performed that confirm this statement, and these results are shown in Table 1.

vMPS | exact | |||||
---|---|---|---|---|---|---|

20 | 12 | 50 | -1.2259 | -1.0804 | -1.0488 | * |

100 | -1.2177 | -1.0646 | -1.03116 | * | ||

* | * | * | -1.0008 | |||

16 | 50 | -0.7972 | -0.5458 | -0.5205 | * | |

100 | -0.7860 | -0.5179 | -0.49513 | * | ||

* | * | * | -0.4639 |

exact | ||||
---|---|---|---|---|

50 | 30 | 50 | -1.2272 | * |

100 | -1.2191 | * | ||

* | -1.0008 | |||

40 | 50 | -0.7974 | * | |

100 | -0.7862 | * | ||

* | -0.4671 |

### iv.2 Correlation functions

Two-particle correlation functions are important quantities in the analysis of lattice systems, because they usually display the physics (for instance the appearance of magnetism) present in the system . In this Section we show that in our approach, these correlation functions are easily extracted from the 2DM, and compare our results to those in Sorella et al. (1990); Ogata and Shiba (1990).

#### Charge correlation

The two-particle charge correlation function is defined as:

(41) |

in which the notation denotes the expectation value. The function is independent of the specific choice of the index because of the periodic boundary conditions. The expression in Eq. (41) can be written in terms of the matrix:

(42) |

and in fact only the singlet part of the matrix appears:

(43) |

In translationally invariant systems one usually takes the Fourier transform of the correlation function,

(44) |

In Figs. 5, 6 and 7, has been plotted for , and half filling respectively, using both and conditions. Comparing the with the results the same trends can be noticed as for the energy and the momentum distributions. For half filling (Fig. 7) the and results are in nice agreement. Moving away from half-filling (Figs. 5 and 6) there is only agreement for small values of . For larger values of strange oscillations appear in the results. So in this limit not only the energy, but the entire physical content of the -2DM cannot be trusted. This is once again an indication that the conditions fail to describe the strong-correlation limit away from half-filling. The results compare well, both in shape and magnitude, with the results from Quantum Monte Carlo Sorella et al. (1990), and the Bethe-ansatz results in the strong-correlation limit Ogata and Shiba (1990).

#### Spin correlation

The two-particle spin-correlation function defined as:

(45) |

can be expressed as a function of the matrix:

(46) |

Written in terms of the spin-coupled matrix, only the triplet part contributes:

(47) |

Fourier transforming yields the momentum dependent spin-correlation function:

(48) |

We have plotted this object in Figs. 8, 9 and 10 for , and half filling respectively, using both as conditions. Unsurprisingly, a good agreement is observed between the and results for small values of . At larger values of , away from half filling, results are very poor with the conditions, especially in the -filled case, where the correlation function is wildly oscillating. More surprising, is that the spin-correlation function for the half-filled lattice in the large- limit is also incorrect in the approximation. In the strong-correlation limit, for half-filling, the spin part of the Hubbard model is identical to the Heisenberg model Ogata and Shiba (1990), for which the spin-correlation function has a singularity at two times the Fermi momentum , which is exactly what we see in the results. Below half-filling the singularity in the large- limit splits and shifts to smaller values of , as obserbed in the figures 8 and 9. This is in agreement with the results in Ogata and Shiba (1990) and Sorella et al. (1990). In conclusion we can say that the 2DM obtained with the conditions, correctly describes the physics that governs the spin-correlation function, whereas the conditions do not. It is also important to note, that even though the results for the energy are good for the half-filled lattice, the 2DM is flawed, because the spin-correlation function is not correctly described.

## V Conclusion

In this article we have studied the one-dimensional Hubbard model at various fillings using the v2DM method with both two- and three-index constraints. We have shown that it is possible to obtain a huge reduction in the computational cost of a basic matrix computation by exploiting all available symmetries, i.e. spin, translation invariance and space-inversion parity. Using this reduction it was possible to compare the computational scaling of different semidefinite programming algorithms with increasing lattice size. We found that, for this particular type of problem, the boundary point method outperforms interior point methods by several orders of magnitude. To gauge the quality of the variationally obtained 2DM we compared several ground-state properties to reference results. We found that, for half filling, the ground-state energy is well described by the two-index conditions. When moving away from half filling, however, we see that the three-index conditions are needed to obtain decent results. An explanation of why this happens has been the subject of a different article Verstichel et al. (2012). The need for three-index constraints was even more obvious when we looked at the spin and charge correlation functions. It was also seen that, even though the energy was relatively correct for the half-filled lattice, the 2DM was flawed, because the spin correlation function was incorrect. This study shows that the exploitation of symmetry opens the possibility for a study of the two-dimensional Hubbard model for relevant lattice sizes. To obtain a decent accuracy, however, it will be necessary to include the three-index constraints, which is computationally hard. One way around this was set forward in Verstichel et al. (2012) with the use of lifting conditions Mazziotti (2002); Hammond and Mazziotti (2005).

## Vi Acknowledgements

We gratefully acknowledge financial support from FWO-Flanders and the research council of Ghent University. B.V., H.V.A., P.B., S.W. and D.V.N. are Members of the QCMM alliance Ghent-Brussels.

## References

- Dirac (1930) P. Dirac, Mathematical Proceedings of the Cambridge Philosophical Society 26, 376 (1930).
- Husimi (1940) K. Husimi, Proc. Phys.-Math. Soc. Japan 22, 264 (1940).
- Fetter and Walecka (1971) L. A. Fetter and J. D. Walecka, Quantum theory of many-particle systems (McGraw-Hill book company, 1971).
- Dickhoff and Van Neck (2008) W. H. Dickhoff and D. Van Neck, Many-Body Theory Exposed! (World Scientific, 2008).
- Löwdin (1955) P. Löwdin, Phys. Rev. 97, 1474 (1955).
- Coleman (2000) A. J. Coleman, Many-electron densities and reduced density matrices, mathematical and computational chemistry (Kluwer Academic/Plenum Publishers, 2000).
- Mayer (1955) J. Mayer, Phys. Rev. 100, 1579 (1955).
- Tredgold (1957) R. H. Tredgold, Phys. Rev. 105, 1421 (1957).
- Coleman (1963) A. J. Coleman, Rev. Mod. Phys. 35, 668 (1963).
- Garrod and Percus (1964) C. Garrod and J. K. Percus, J. Math. Phys. 5, 1756 (1964).
- Garrod and Fusco (1976) C. Garrod and M. A. Fusco, Int. J. Quantum Chem. 10, 495 (1976).
- Garrod et al. (1975) C. Garrod, M. V. Mihailovic, and M. Rosina, J. Math. Phys. 16, 868 (1975), URL http://link.aip.org/link/?JMP/16/868/1.
- Mihailovic and Rosina (1974) M. Mihailovic and M. Rosina, Nuclear Physics A 237, 221 (1974).
- Rosina and Garrod (1975) M. Rosina and C. Garrod, J. Comp. Phys. 18, 300 (1975).
- Nakata et al. (2001) M. Nakata, H. Nakatsuji, M. Ehara, M. Fukuda, K. Nakata, and K. Fujisawa, J. Chem. Phys. 114, 8282 (2001).
- Mazziotti (2002) D. A. Mazziotti, Phys. Rev. A 65, 062511 (2002).
- Zhao et al. (2004) Z. Zhao, B. J. Braams, M. Fukuda, M. L. Overton, and J. K. Percus, J. Chem. Phys. 120, 2095 (2004).
- Hammond and Mazziotti (2005) J. R. Hammond and D. A. Mazziotti, Phys. Rev. A 71, 062503 (2005).
- Nakata et al. (2008) M. Nakata, B. J. Braams, K. Fujisawa, M. Fukuda, J. K. Percus, M. Yamashita, and Z. Zhao, J. Chem. Phys. 128, 164113 (2008).
- Mazziotti (2005) D. A. Mazziotti, Phys. Rev. A 72, 032510 (2005).
- Gidofalvi and Mazziotti (2006) G. Gidofalvi and D. A. Mazziotti, J. Chem. Phys. 125, 114102 (2006).
- Mazziotti (2007) D. A. Mazziotti, Reduced-Density-Matrix Mechanics: With Aplication to Many-Electron Atoms and Molecules, vol. 134 (Wiley: New York, 2007).
- Braams et al. (2007) B. J. Braams, J. K. Percus, and Z. Zhao, Reduced-Density-Matrix Mechanics: With Aplication to Many-Electron Atoms and Molecules, vol. 134 (Wiley: New York, 2007).
- Van Neck and Ayers (2007) D. Van Neck and P. W. Ayers, Phys. Rev. A 75, 032502 (2007).
- Verstichel et al. (2010) B. Verstichel, H. van Aggelen, D. V. Neck, P. W. Ayers, and P. Bultinck, J. Chem. Phys. 132, 114113 (2010).
- Mazziotti (2004) D. A. Mazziotti, Phys. Rev. Lett. 93, 213001 (2004).
- Mazziotti (2011) D. A. Mazziotti, Phys. Rev. Lett. 106, 083001 (2011).
- Verstichel et al. (2011) B. Verstichel, H. van Aggelen, D. V. Neck, P. Bultinck, and S. D. Baerdemacker, Comput. Phys. Commun. 182, 1235 (2011).
- Hubbard (1963) J. Hubbard, Proc. R. Soc. Lond. A 276, 238 (1963).
- Sorella et al. (1990) S. Sorella, A. Parola, M. Parrinello, and E. Tosatti, EPL (Europhysics Letters) 12, 721 (1990).
- Ogata and Shiba (1990) M. Ogata and H. Shiba, Phys. Rev. B 41, 2326 (1990), URL http://link.aps.org/doi/10.1103/PhysRevB.41.2326.
- Verstichel et al. (2009) B. Verstichel, H. van Aggelen, D. Van Neck, P. W. Ayers, and P. Bultinck, Phys. Rev. A 80, 032508 (2009).
- Hammond and Mazziotti (2006) J. R. Hammond and D. A. Mazziotti, Phys. Rev. A 73, 062505 (2006).
- Shenvi and Izmaylov (2010) N. Shenvi and A. F. Izmaylov, Phys. Rev. Lett. 105, 213003 (2010).
- Verstichel et al. (2012) B. Verstichel, H. van Aggelen, W. Poelmans, and D. Van Neck, Phys. Rev. Lett. 108, 213001 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.108.213001.
- Bethe (1931) H. Bethe, Zeitschrift für Physik A Hadrons and Nuclei 71, 205 (1931), ISSN 0939-7922, URL http://dx.doi.org/10.1007/BF01341708.
- Lieb and Wu (1968) E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968), URL http://link.aps.org/doi/10.1103/PhysRevLett.20.1445.
- Essler et al. (1991) F. H. L. Essler, V. E. Korepin, and K. Schoutens, Phys. Rev. Lett. 67, 3848 (1991), URL http://link.aps.org/doi/10.1103/PhysRevLett.67.3848.
- Essler et al. (2005) F. H. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin, The One-Dimensional Hubbard Model (Cambridge University Press, 2005).
- Schollwöck (2011) U. Schollwöck, Annals of Physics 326, 96 (2011), ISSN 0003-4916, URL http://www.sciencedirect.com/science/article/pii/S0003491610001752.
- Verstraete et al. (2008) F. Verstraete, V. Murg, and J. Cirac, Advances in Physics 57, 143 (2008).
- Chan and Zgid (2009) G. K.-L. Chan and D. Zgid (Elsevier, 2009), vol. 5 of Annual Reports in Computational Chemistry, pp. 149 – 162.
- Wouters et al. (2012) S. Wouters, P. A. Limacher, D. Van Neck, and P. W. Ayers, J. Chem. Phys. 136, 134110 (2012), eprint cond-mat.str-el/1202.0177v2.
- Lieb and Wu (2003) E. H. Lieb and F. Wu, Physica A: Statistical Mechanics and its Applications 321, 1 (2003), ISSN 0378-4371, URL http://www.sciencedirect.com/science/article/pii/S0378437102017855.
- Nakata and Yasuda (2009) M. Nakata and K. Yasuda, Phys. Rev. A 80, 042109 (2009), URL http://link.aps.org/doi/10.1103/PhysRevA.80.042109.
- van Aggelen et al. (2009) H. van Aggelen, P. Bultinck, B. Verstichel, D. Van Neck, and P. W. Ayers, Phys. Chem. Chem. Phys. 11, 5558 (2009).
- Krivnov et al. (1990) V. Krivnov, A. Ovchinnikov, and V. Cheranovskii, Teoreticheskaya i Matematicheskaya Fizika 82, 216 (1990).