Correlation density matrices for 1- dimensional quantum chains based on the density matrix renormalization group

A useful concept for finding numerically the dominant correlations of a given ground state in an interacting quantum lattice system in an unbiased way is the correlation density matrix. For two disjoint, separated clusters, it is defined to be the density matrix of their union minus the direct product of their individual density matrices and contains all correlations between the two clusters. We show how to extract from the correlation density matrix a general overview of the correlations as well as detailed information on the operators carrying long-range correlations and the spatial dependence of their correlation functions. To determine the correlation density matrix, we calculate the ground state for a class of spinless extended Hubbard models using the density matrix renormalization group. This numerical method is based on matrix product states for which the correlation density matrix can be obtained straightforwardly. In an appendix, we give a detailed tutorial introduction to our variational matrix product state approach for ground state calculations for 1- dimensional quantum chain models. We show in detail how matrix product states overcome the problem of large Hilbert space dimensions in these models and describe all techniques which are needed for handling them in practice.


Introduction
In an interacting quantum lattice model the ground state may have several kinds of correlations, such as long-range order, power-law, or exponentially decaying correlations. In the numerical treatment of such a model it is not clear a priori what kind of correlation will be dominant and what kind of operators corresponds to these correlations. Before calculating correlation functions, one typically chooses in advance which operators to consider, using prior knowledge and making initial assumptions. The need to make such choices introduces a certain bias into the investigation, which can be somewhat unsatisfying, especially when hidden or exotic correlations are present.
which is completely unbiased except for the specification of the clusters. If the two clusters were not correlated at all, this would implyρ AB =ρ A ⊗ρ B and thereforê ρ C = 0. The CDM encodes all possible correlations between the clusters A and B, as can be seen from the fact that whereÔ A andÔ B are operators acting on clusters A and B, respectively.

Lessons from Luttinger liquid theory
To extract useful information from the CDM, it will be helpful to develop some intuition for its general structure. To this end, let us recall some fundamental facts from onedimensional critical fermion systems. They are described by Luttinger liquid theory, in which one of the key parameters is the Fermi wave vector k F . The asymptotic behavior of any kind of correlation or Green's function is typically an oscillation inside a power-law envelope, for some exponent γ, where m is some integer. For the particular model to be used in this study, a nontrivial mapping is known to a free fermion chain [2], a special case of Luttinger liquid. Renormalization group theory [6] quite generally implies the existence of scaling operators in any critical system such as a Luttinger liquid. They are eigenvectors of the renormalization transformation and consequently their correlations are purely of a form like (1.3) for all r, not just asymptotically. The scaling operators usually have complicated forms. The correlation of a simple operator (e.g. fermion density n(x) at position x along a chain) has overlap with various scaling operators, and correspondingly the correlation function of that simple operator is a linear combination of contributions like (1.3) from those scaling operators.
Our aim is to discover the leading scaling operators numerically. The leading scaling operator encodes all the local fluctuations that are correlated with faraway parts of the system. Intuitively, for a given cluster A, that operator does not depend significantly on the exact position of the (distant) cluster B. That is particularly obvious in a one dimensional system: any correlation at distances r > r must be propagated through some sort of correlation at r, so we expect the same operators from cluster A to be involved inρ C (r), irrespective of the distance r.
This suggests an ansatz for leading contributions in the CDM: HereÔ A,s andÔ B,s are a pair of (distance-independent) scaling operators acting on clusters A and B, respectively, k s is the characteristic wave vector for oscillations in their correlation, and γ s is the corresponding scaling exponent. When k s = 0, the operator pairs must themselves come in pairs, labelled, say, by s and s + 1, with k s+1 = −k s , c s+1 = c * s , and γ s+1 = γ s , so thatρ C is hermitian. The scaling operators for each cluster form an orthonormal set. We expect that only a few terms in the sum in (1.4) capture most of the weight. Correspondingly, it may be feasible to truncate the complete basis setsÔ A,s andÔ B,s to a smaller set of "dominant operators", whose correlators carry the dominant correlations of the system. The ansatz (1.4) will guide our steps in the ensuing analysis; at the end, we shall check how well it is satisfied by the actual CDMs calculated for the model studied in this paper (see section 6.1.2).
Notice that although a particular correlation function may have nodes, see (1.3), for a CDM of the form (1.4) the norm, is monotonically decaying with r. This expresses the fact that information can only be lost with increasing distance, never restored, in a one-dimensional system.

Operator basis and f-matrix
In [1] the operators entering the dominant correlation were found by a kind of singular value decomposition (SVD), which was done independently for each separation. However, the operators obtained from the SVD will in general be different for different separations r. This does not correspond to the form (1.4), where the operators are distance-independent and only the coefficients are r-dependent. Therefore, we shall explore in this paper a new scheme to decompose the CDMs for all separations in concert, so as to obtain a small set of scaling operators characterizing the dominant correlations at any (sufficiently large) separation. We decomposeρ C in the form where the S i represent the symmetry-sectors of the discrete, Abelian symmetries of the Hamiltonian (see section 3.3). The subscript of the brackets indicates that the decomposition within the brackets is done for each symmetry-sector individually. This decomposition is possible for any two complete, r-independent operator setsÔ A,µ andÔ B,µ acting on the part of the Hilbert space of clusters A and B, respectively, which correspond to the symmetry sector S i . The goal is to find two operator setŝ O A,µ andÔ B,µ such that these operator sets may be truncated to a small number of operators each, while still bearing the dominant correlations of the system. The distance dependence of the CDM is then only contained in the matrix f µ,µ (r). Then, all analysis concerning the distance-dependence of correlations can be done in terms of this f-matrix.

Ground state calculation with DMRG
The CDM in [1] was calculated using the full ground state obtained from exact diagonalization. This limits the system size, so that the method was appropriate mainly in cases of rapidly decaying, or non-decaying correlations -not for critical or slowly decaying ones. In the present work, we use the density matrix renormalization group (DMRG) [3] (see the excellent review by U. Schollwöck [4]) to compute the ground state for a ladder system which is known to have algebraic correlations [2]. We use the matrix product state (MPS) formulation of DMRG [5] in which an efficient variational procedure is used to obtain the ground state.

Structure of the paper
The structure of the main body of the paper is as follows: in section 2 we introduce the model to be considered for explicit calculations. In section 3 we show how the CDM is defined, how to calculate it, and explain how a first overview of the relative strengths of various types of correlations can be obtained. In section 4 we show how to analyze the CDM and its distance dependence. Sections 5 to 7 present our numerical results, and section 8 our conclusions. In an extended appendix we offer a tutorial introduction to the MPS formulation of DMRG, and also explain how it can be used to efficiently calculate the CDM.

Model
To be concrete in the following analysis of the CDM, we begin by introducing the model for which we did our numerical calculations. This model contains rich physics and its treatment below can readily be generalized to other models.

Definition of the model
We analyze the CDM for a class of spinless extended Hubbard models for fermions, which was intensely studied by Cheong and Henley [2]. They computed correlation functions up to separations of about r = 20, using nontrivial mappings to free fermions and hardcore bosons. The correlation functions are calculated with an intervening-particle expansion [2], which expresses the correlation functions in terms of one-dimensional Fermi-sea expectation values (an evaluation of the CDM for that model has also been done by Cheong and Henley [1], using exact diagonalization, but the system sizes are too short to be conclusive). For spinless fermions on a two-leg ladder with length N , we use the following Hamiltonian:  whereĉ a,x destroys a spinless fermion on leg a and rung x, andn a,x =ĉ † a,xĉ a,x is the corresponding number operator. Effectively, the model corresponds to a one-dimensional pseudo-spin chain, where the a = 1 leg is denoted by spin ↑ and the a = 2 leg by spin ↓. Hence, in the following sections which generally apply to quantum chain models we will treat this model as a quantum chain consisting of N sites and return to view the system as a ladder model in the sections where we discuss our results.
We will focus on infinite nearest-neighbour repulsion V → ∞, which we treat differently along the legs and the rungs in our numerical calculations. In the pseudospin description we can enforce the nearest-neighbour exclusion along rungs by removing double occupancy from the local Hilbert space of the pseudo-spin sites. The nearestneighbour exclusion along the legs cannot be implemented so easily and we mimic V → ∞ by a value of V which is much larger than all the other energies in the Hamiltonian (typically V /t = 10 4 ).
For fermionic systems, the fermionic sign due to the anti-commutation relations of the fermionic creation-and annihilation-operators needs to be taken into account. Specifically, we have to choose an order in which we pick the Fock basis, where we have to keep in mind that this choice produces a so called Jordan-Wigner-string of the form x −1 x =x+1 e iπn x when evaluating correlators ĉ xĉ † x at distance r = |x − x |. In the present system it is convenient to choose this order such that the operators of the two sites of a rung are succeeding each other (see figure 1), as this choice yields the shortest Jordan-Wigner strings.

Expectations for simple limiting cases
Setting t ≡ 1 as a reference scale, we are left with two parameters in the Hamiltonian: the rung hopping t ⊥ and the correlated hopping t c . The physics of the system is governed by the competition of t ⊥ to localize the fermions on the rungs and t c to pair the fermions. There are three limiting cases which have been studied in detail by Cheong and Henley [1,2].
(i) The paired limit, t c t , t ⊥ (we used t c /t = 10 2 and t ⊥ = 0 for our calculations).
In this limit the fermions form tight pairs which behave similar to hardcore bosons [2]. For two given rungs x and x + 1, there are two possibilities to create a pair of fermions, due to infinite nearest-neighbour repulsion:ĉ † ↑xĉ † ↓x+1 andĉ † ↓xĉ † ↑x+1 . It has been shown in [2] that, based on these two bound pairs, one may classify the bound pairs in two flavours along the ladder and that the ground state has only one definite flavour, causing a twofold symmetry breaking in the ground state. This symmetry breaking introduces complications that will be addressed below. The dominant correlations are expected to be charge-density correlations at short distances and two-particle at long distances. These charge-density and two-particle correlations decay as power laws, oscillating with k = 2k F , where the Fermi wavelength k F is related to the filling as k F = 2ν [2]. In this system, the one-particle correlations are suppressed and are expected to decay exponentially, as a nonzero expectation value depends on a local fluctuation completely filling the rungs between the clusters (as elaborated in section 6.2).
(ii) The two-leg limit, t ⊥ t , t c = 0. In this limit the two legs are decoupled with respect to hopping, but still the infinite nearest-neighbour repulsion introduces correlations between the two legs. At large distances, power-law charge-density correlations dominate, while two-particle correlations show much faster power-law decay and one-particle correlations decay exponentially.
(iii) The rung-fermion limit, t ⊥ t , t c = 0. In this limit the particles are delocalized along the rungs. For fillings smaller than quarter-filling, charge-density , oneparticle and two-particle correlations all decay as power laws where charge-density correlations dominate at large distances.
Our analysis in this paper is limited to the case (i), where DMRG also showed best performance.

Smooth boundary conditions
For a ladder of length N (treated as a pseudo-spin chain), we have attempted to reduce effects from the boundaries by implementing smooth boundary conditions, adopting a strategy proposed in [7] for a spin chain to our present fermionic system. (Alternatively, it is possible to use periodic boundary conditions [5]. However, this leads to some difficulties, since it is not possible to work with orthonormal basis sets describing the left or right part of the chain with respect to a given site.) Smooth boundary conditions are open boundary conditions together with an artificial decay of all terms of the Hamiltonian over the last M rungs at each end of the chain. We shall calculate expectation values only of operators located in the central part of the system (sites x, with M < x ≤ N − M ), thus the system's effective length is N = N − 2M .
For both smooth and open boundary conditions the average site filling strongly decreases near the boundaries. To determine the average filling ν, which influences the system's correlations in an important manner, we thus use only the central N sites: Due to the infinite nearest neighbour repulsion, this implies that ν ∈ [0, 0.5].

Calculation of the CDM
Throughout the paper we will use the Frobenius inner product and norm for any matrices M ij and M ij of matching dimension,

Definition of the CDM
We take two disjoint, separated clusters A and B of equal size from a one-dimensional quantum chain, i.e. two sets of adjacent sites x A 1 , . . . , x A n and x B 1 , . . . , x B n where n is the size of the clusters and all the indices x are distinct from each other. The local Hilbert spaces of clusters A and B with dimension d n are described in terms of sets of basis states |α and |β , which are product states of the local states of each site in the cluster. The CDM of the two clusters, defined by (1.1), can be expanded in this basis aŝ For processing the CDM we fuse the two indices of each cluster [1]: withα = (αα ) andβ = (ββ ), and denote the reshaped objectρ C itself by an extra tilde. This corresponds to a partial transpose of the CDM (note thatρ C is no longer a symmetric tensor). For the CDM expressed in the indicesα andβ, we may use the Frobenius inner product (3.1) and norm (3.2).
To study the distance dependence of the correlations, we vary the position of the clusters A and B, resulting in a position-dependent CDMρ C x A 1 , x B 1 . If the system is translationally invariant, this object depends only on the distance r = |x A 1 − x B 1 | (the minimal distance for two adjacent clusters is equal to the cluster size n). For a finite system, though,ρ C will also depend on 1 2 x A 1 + x B 1 , at best weakly if the system is long. Strategies for minimizing the dependence on 1 2 by taking suitable averages will be discussed in section 3.4.

DMRG-calculation of the CDM
The fact that the Hamiltonian in (2.1) is a one-dimensional pseudo-spin chain allows us to calculate ground state properties with the density matrix renormalization group (DMRG) [3]. Using the variational matrix product state formulation of that method (see appendix for a detailed description), we calculated the ground state of the Hamiltonian in (2.1) for several values of t ⊥ and t c . The framework of MPS also allows the CDM to be calculated efficiently (see section A.2.7 for details). Limiting ourselves to the case t ⊥ = 0 in this paper, we have calculated the CDM derived from the ground state for distances up to 40 rungs, which is significantly larger than in previous approaches.

Symmetry sectors
All the symmetries of the Hamiltonian are reflected in the CDM, making the CDM block-diagonal, where each block can be labeled uniquely by a set of quantum numbers that are conserved by the Hamiltonian. This means for Abelian symmetries (which are the only ones we are considering in practice), that the CDM in the original form ρ C αβ,α β fulfills Q α + Q β = Q α + Q β , where Q α corresponds to the quantum numbers of state |α , etc. The rearrangement of the CDM intoρ C αβ then implies ∆Qα = −∆Qβ with ∆Qα ≡ Q α − Q α and ∆Qβ ≡ Q β − Q β . Sinceρ AB is hermitian, for every block of the CDM involving ∆Qα (∆Qβ) there has to be a block involving −∆Qα (−∆Qβ), respectively. Therefore, it is convenient to sort the various parts of the CDM in terms of their change in quantum numbers ∆Q ≡ |∆Qα| = |∆Qβ| and to analyze each symmetry sector individually.
To obtain a general classification of the CDM we sort the various contributions of the CDM according to the conserved quantum number(s) Q. In the case of the Hamiltonian in (2.1), we consider particle conservation (Q =N tot ) which breaks the CDM into blocks with well-defined particle transfer ∆N ≡ |∆Nα| = |∆Nβ| between clusters A and B. The following r.m.s. net correlations then is a measure for the correlations with transfer of ∆N particles between A and B (with ∆N = 0, 1, 2): Here the notationα ≡ (αα ) ∈ S ∆N indicates that only pairs of states (αα ) are considered which differ by ∆N in particle number (similarly forβ ≡ (ββ ) ∈ S ∆N ). In the following we will call correlations involving ∆N = 0, 1, 2 particles charge-density correlations (CD), one-particle correlations (1P), and two-particle correlations (2P), respectively. The following analysis is done for each symmetry sector individually. Depending on the decay of the r.m.s. net correlations (3.5), some symmetry sectors may become irrelevant with increasing distance.

"Restoration" of numerically broken symmetries
Although we have tried to minimize the effect of boundaries, our numerical methods for calculating the ground state and CDM do not produce strictly translationally invariant results. (In contrast, analyses based on exact diagonalization start from a ground state wavefunction in which the symmetry (in a finite system) is restored, even if there is a symmetry breaking in the thermodynamic limit.) Therefore, we construct the CDM ρ C (r) for a given distance r from an average over several CDMsρ C (x, x ) with constant r = |x − x |, where x and x give the position of the first site of clusters A and B, respectively.
Moreover, if the exact ground state is degenerate under a discrete symmetry, we expect that DMRG breaks this symmetry unless it is implemented explicitly in the code. As mentioned in section 2.2 for the specific models of this paper we expect a discrete symmetry under interchange of legs for some parameter regimes. Since we did not implement this symmetry explicitly in our code, we also average the CDM by interchanging the legs of the ladder. Thus, all the data analysis presented in subsequent sections will be based on using the following "symmetry-restored" form of the CDM, whereρ C is obtained fromρ C by interchanging the legs of the ladder, and N is some normalization factor. One might argue that it is not sufficient to average over the broken symmetry w.r.t. leg-interchange on the level of the density matrix, but that instead the symmetry should be restored on the level of the ground state wave function. Specifically, for a ground state |ψ 1 (however it is calculated) which breaks this symmetry, we could restore the symmetry in the following way, where |ψ 2 =Ŝ |ψ 1 andŜ describes the action of interchanging the legs. This would lead to a total density matrix Now, for two clusters A and B, the first two terms on the r.h.s. yield the CDM of (3.6), while the last two terms turn out to be negligible when traced out over all sites except for the two local clusters A and B. This follows from |ψ 1 and |ψ 2 being orthogonal, hence tr(|ψ 1 ψ 2 |) = ψ 2 |ψ 1 = 0, implying that for a long chain with local clusters A and B, the reduced density matrixρ AB,12 ≡ tr x / ∈A,B (|ψ 1 ψ 2 |) will be very close to zero due to the orthogonality of the wave functions on the sites outside of clusters A and B. Consequently, it is sufficient to retain only the first two terms of (3.8), i.e. to restore the broken symmetry on the level of the density matrices only, as done in (3.6).

Finding a distance-independent operator basis
The goal of this section is to extract a (likely) small set of operators from the CDM, which will describe the dominant correlations in the system as a function of distance. We will assume in this section that the CDM does not include any broken symmetries as indicated in section 3.4.

Need for operator bases for clusters A and B
As already mentioned, the CDM (obtained from (3.6)) may be investigated by applying a singular value decomposition (SVD) for each distance individually [1]: 1) or, in operator notation: = O B,s ββ form a complete set in the operator space of clusters A and B, respectively, using the inner product as in (3.1). The set includes operators with w s = 0, such as the identity operator, since these will be produced by the SVD. The SVD (4.2) yields for each specific distance r a set of operatorsÔ A,s (r) andÔ B,s (r) acting on clusters A and B, respectively.
However, the dominant operators so obtained, i.e. the ones with large weight from the SVD ofρ C (r), are likely not the same as each other for different distances and hence not convenient for characterizing the "dominant correlations" of the system. What is needed, evidently, is a strategy for reducing the numerous sets of operatorsÔ A,s (r) and O B,s (r) to two "basis sets of operators" for clusters A and B, respectively, sayÔ A,µ and O B,µ , which are r-independent and whose correlators yield the dominant correlations in the system in the spirit of (1.4). (For a translationally invariant system the two sets have to be equal for both clusters A and B, but we will treat them independently in the analysis.) Following the ansatz (1.4) from the Luttinger liquid theory, these operators ought to be distance-independent, carrying common correlation content for all distances. Thus we seek an expansion ofρ C (r) of the form (1.6), in which only the coefficients, not the operators, are r-dependent.

Construction of operator bases
We have explored a number of different strategies for extracting operators from the CDM which carry common information for all distances. We will discuss in detail only one of these, which is rather simple to formulate and reliably yields operator sets with the desired properties. (Several other strategies yielded equivalent results, but in a somewhat more cumbersome fashion.) The simplest possible strategy one may try is to average over all the CDMs at different distances and to singular-value decompose the resulting crude "average CDM". However, since the elements for the CDM are expected to be oscillating functions of r, such a crude average can cancel out important contributions of the CDM. Thus we need a procedure that avoids such possible cancellations. To this end, we construct the following operators, bilinear in the CDM: with matrix elements We normalize by ρ C (r) 2 in order to treat the operator correlations ofρ C (r) for different distances on an equal footing. Note that the eigenvalue decomposition on the hermitian matrices K A (r) and K B (r) (in short K-matrices) yields the same operatorŝ In particular, it no longer contains any oscillating parts (in contrast to (1.4)), and hence is suitable for being averaged over r.
Summing up the K X -matrices over a range R of distances (r ∈ R, where R will be specified below) gives a meanK X -matrix for cluster X (= A, B), namelȳ K X,R ≡ r∈RK X (r). We do not divide the latter expression by the number of terms in the sum (as would be required for a proper mean), as at this stage we are only interested in the operator eigendecomposition, with the operators normalized such that Ô X,R,µ = 1. The operator setÔ X,R,µ gives an orthonormal, r-independent basis for cluster X. In practice, however, many of the w R,µ (which turn out to be the same for X = A or B) will be very small. Thus, it will be sufficient to work with a truncated set of these operators having significant weight. To explore the extent to whichK X depends on the summation range, we shall study several such ranges: R all includes all distances, R short short distances (first third of distances analyzed), R int intermediate distances (second third) and R long long distances (last third). The resulting (truncated) sets of operators can be compared via their mutual overlap matrix O RR µµ = tr(Ô R,X,µÔR ,X,µ ), or more simply, by the single number O RR = µµ (O RR µµ ) 2 , which may be interpreted as the dimension of the common subspace of the two operator sets. The value of O RR ranges from 0 to dim(Ô R,X,µ ). By comparing O RR for the different distance ranges, additional clues can be obtained about how the relative weight of correlations evolves from short to long distances. (Such a comparison is carried out in table 1 below.)

Definition of f-Matrix
Once a convenient basis of operatorsÔ A,µ andÔ B,µ has been found, the correlation density matrix can be expanded in terms of this basis as in (1.6), with matrix elements For complete operator spacesÔ A,µ andÔ B,µ , by definition, the set of amplitudes squared sum up to the norm of the CDM: However, as alluded to above, we expect that the dominant correlators can be expressed in terms of a truncated set of dominant operators. If the sum on the left hand side of (4.9) is restricted to this truncated set, its deviation from the right hand side gives an estimate of how wellρ C is represented by the truncated set of operators. It will turn out that only a handful of dominant operators (typically 4 or 6) are needed, implying very significant simplifications in the analysis. Thus, the data analysis will be done in terms of the matrices f µ,µ (r) (in short "f-matrix") for this truncated set of dominant operators.

Fourier-analysis and decay of f-matrix
According to the expectations expressed in (1.4), the elements of the f-matrix are expected to be products of oscillating and decaying functions of r. The corresponding dominant wave vectors can be identified via Fourier transform on each element of the f-matrix. For an oscillating function times a monotonically decaying envelope, the peaks of the Fourier spectrum of the oscillating function will be broadened by the presence of the envelope. To minimize this unwanted broadening, we introduce a rescaled fmatrix (denoted by a tilde),f µ,µ (r) = u (r) f µ,µ (r), where the positive weightingfunction u (r) is chosen such that all values of |f µ,µ (r) | are of the same order, and Fourier decompose the rescaledf -matrix asf µ,µ (k) = r e −ikrf µ,µ (r). Its norm f (k) 2 = µµ |f µ,µ (r) | 2 , plotted as a function of k, will contain distinct peaks that indicate which wave vectors characterize the dominant correlations. Subsequently, the elements of the f-matrix, can be fitted to the forms µ,µ are complex amplitudes, f j (r) describes the decay with distance (e.g. f j (r) = r −γ j or e −r/r j for power-law or exponential decay, respectively), and k j is a set of dominant wave vectors. The latter appear pairwise in combinations (+k; −k), The results of such a fit for each pair of dominant operatorsÔ A,µ andÔ B,µ , is the final outcome of our analysis, since it contains the information needed to check the applicability of ansatz (1.3).

Numerical results: general remarks
In this section, we illustrate the analysis proposed above for the model introduced in section 2. We will focus on the limiting case of large t c , which we expect to have the most complex behavior among all three limiting cases introduced in [1] and [2]. After some preliminary analysis, we will discuss in section 6 each of the three symmetry sectors (CD, 1P, and 2P) characterized by the operators' fermion number, and in section 7 compare our results to those found by [2] using a different method.

Specification of the clusters A and B
For the following analysis it is convenient to take the size of the clusters A and B to be two rungs, because clusters of at least that size allow for up to two particles in one cluster (due to infinite nearest-neighbour repulsion). Thus, correlations involving ∆N = 0, 1, 2 are possible, i.e CD, 1P, and 2P correlations, respectively. Note that larger clusters can be studied, but would significantly increase numerical costs. Taking into account the infinite nearest-neighbour repulsion, clusters of size two have a seven-dimensional Hilbert space spanned by the kets |00 , |0 ↑ , |0 ↓ , |↑ 0 , |↓ 0 , |↑↓ , |↓↑ , where the first (second) entry corresponds to the first (second) rung, 0 represents an empty rung and ↑ and ↓ a fermion on the upper and lower leg in pseudo-spin notation (recall that we are dealing with spinless fermions). The space of operators acting on a cluster has dimension 7 2 = 49, where the subspaces for ∆N = 0, 1 or 2 have dimensions 21, 24 and 4, respectively, as depicted schematically in figure 2.

Average site occupation
As a first check of the influence of the boundaries, we investigate the average site occupation on the ladder. It is expected to be uniform in a translationally invariant system. However, there are two ways in which our calculation breaks translational symmetry, which cause residual oscillations in the density of particles along the ladders.
Firstly, there is the spontaneous breaking of the pair flavor symmetry described in section 2.2. In the ground state produced by DMRG, all pairs have the same flavor, so only one of the two sublattices actually has any fermions on it. Thus a strong alternation in the density is observed between one leg for even rungs and the other leg for odd rungs; this can be taken care of by the symmetrization with respect to legs (as in (3.6)). Secondly, translational symmetry is broken due to finite size in the DMRG calculation. This induces oscillations in the average occupation as a function of x (see figure 3), whose period is clearly dependent on the filling. In fact, their period is 2k F , so they may be interpreted as Friedel-like oscillations caused by the boundaries. Although the amplitude of density oscillation appears rather flat in the central portion of the system, it does have a minimum there; so we expect that the amplitude in the center of the system would vanish in a sufficiently large system.
Although the intent of the smooth boundary conditions is to minimize effects such as these oscillations, in fact, their amplitude appeared to be of about the same strength independent of whether we used smooth or plain open boundary conditions. We suspect, however, that the amplitude could be reduced by further careful optimization (not attempted here) of the parameters of the smooth boundary conditions.

r.m.s. net correlations w ∆N (r)
The next basic step is to identify the leading correlations in terms of the r.m.s. net correlations w ∆N defined in (3.5). These reveal which sectors of correlations dominate at large distances. The results (see figure 4) show that the r.m.s. net correlations decay exponentially in the 1P sector, whereas they decay algebraically in both the CD and 2P sectors, consistent with [2]. The latter two correlations are comparable in size over a significant range of distances, but for the fillings we investigated, 2P correlations ultimately dominate over CD correlations at the largest distances. Both the CD and 2P r.m.s. net correlations can be fitted to power laws, with the exponent dependent on the filling. The r.m.s. net correlations in each sector are monotonic and only weakly modulated, even though the dominant correlation functions and the dominant parts of the CDM itself are oscillating (as will be discussed in more detail in section 6.1, see, e.g., figure 7). This implies that the correlations in each sector can be represented by a linear combination of correlation functions (associated with different operators) which oscillate out of phase, in such a way that in the sum of their squared moduli the oscillations more or less average out, resulting in an essentially monotonic decay with r, as expected according to (1.5).
We will next apply the analysis proposed in section 4.2 to the respective symmetry sectors (which will provide more exact fits of the exponents of the power-law decays). The analysis in any sector consists of two stages. First, following section 4.2, we try to find an optimal truncated basis which describes best the dominant correlations. Second, we examine the f-matrix of section 4.3 (i.e. represent the CDM in the truncated basis) to see the nature of its r dependence, and to fit this to an appropriate form, following section 4.4.
6. Numerical results: symmetry sectors 6.1. Charge-density correlations 6.1.1. Operator basis First we calculated the mean K-matricesK A,R andK B,R from ρ C R defined in (4.3a) and (4.3b), and obtained operator sets from their eigenvalue The symmetry sectors are ∆N = 0 (blue, no particle transfer, CD), ∆N = 1 (green, transfer of one particle, 1P) and ∆N = 2 (red, transfer of two particles, 2P). We see that CD and 2P correlations decay as power-laws (r −γ , blue and red solid lines) with small residual oscillations at k = 2k F , while the 1P correlations show exponential decay (e −r/r1 , see semi-logarithmic plot in the inset). The value r 1 0.5 for both fillings is reasonable as we would expect a value of the order of one, which is the size of the bound pairs. decomposition, using various distance ranges.
In order to decide how many operators to include in the truncated basis, we used the diagnostic described in section 4.2. In presenting the results, we limit ourselves to cluster A as the results for cluster B are completely analogous. The operator set O A,R all ,µ corresponding to the full range of distances R all (specified in section section 4.2) is used as a reference set to be compared with the operator sets obtained from R short , R int and R long . The results are given in table 1. We see that, for intermediate or long distances, the effective dimension (O R all R int and O R all R long ) of the common operator space shared between the operator setÔ A,R all ,µ and the operator setsÔ A,R int ,µ andÔ A,R long ,µ , respectively, saturates at six even if a larger operator space is allowed. Similarly, also the short-distance operator setÔ A,R short ,µ agrees best with the other three operator sets at dimension six: a further increase of the number of operators, however, adds only operators in the short range sector of the CDM. Hence we truncate to a six-dimensional operator basis. Within this reduced operator space, all dominant correlations are wellcaptured, as can be seen from the relative weights of table 1. For the resulting truncated basis set equation (4.9) holds up to a relative deviation of the order of O (10 −5 ).
Investigating the six-dimensional set of operators in more detail reveals that they can be classified with respect to their symmetry with respect to interchanging the legs of the ladder, i.e. they obeyŜÔ A,R all ,µ = ±Ô A,R all ,µ , withŜ describing the action of interchanging legs. The set breaks into two subsets of three operators each, which have positive or negative parity with respect toŜ, respectively. It turns out that all six operators are linear combinations of operators having matrix elements on the diagonal only, in the representation of figure 2. Moreover, together with the unit matrix they span  the full space of diagonal operators (therefore the dimension of 6 = 7 − 1). Explicitly, the symmetric operators are given bŷ [−6n 0,xn0,x+1 + (n 0,xn↑,x+1 +n ↑,xn0,x+1 +n ↑,xn↓,x+1 + leg symmetrized)](6.1c) and the antisymmetric operators bŷ . We use this operator basis for both cluster A and cluster B.
If we calculate the f-matrix (4.7) based on these operators we see that it breaks into two blocks corresponding to their symmetry with respect to leg interchange.

f-matrix elements: oscillations and decay
We now turn to extracting the distance-dependence of the dominant correlation in this symmetry sector, which is now visualizable since we drastically reduced the operator space to six dimensions. All relevant information is contained in the f-matrix and its Fourier transform. The first step is to identify the oscillation wave vector(s) k to be used as initial guesses in the fit. A general method is to plot the Fourier spectrum f (k) of the rescaled fmatrix ( figure 5). When using a logarithmic scale for the vertical axis, even sub-leading contributions show up clearly. We find that the spectra belonging to the symmetric and anti-symmetric operators are shifted against each other by π. This relative phase Figure 5. Fourier transform of the rescaled f-matrixf for CD correlations based on operators chosen from a reduced six-dimensional operator space, for a filling of (a) ν = 0.248 and (b) ν = 0.286. We obtain these Fourier spectra from the rescaled f-matrixf µ,µ (r) = r γ f µ,µ (r), with γ extracted from a power-law fit on |f µ,µ (r) |. The Fourier spectrum breaks up into a contribution coming from the operators symmetric or antisymmetric under leg-interchange, labelledf + (blue) and f − (red), respectively. The spectrum off + shows strong peaks at k = ±2k F (dashed lines) and a smaller peak at k = 0 with k F /π = ν. The spectrum off − , having peaks at k = ±2k F + π (dashed lines) and k = π, is shifted w.r.t.f + by π. For a filling close to 1 4 the dominant peaks off ± , at k = ±2k F and k = ±2k F + π. are nearly at the same position.
shift implies a trivial additional distance dependence of e iπr of f − (r) with respect to f + (r), reflecting the different parity under leg interchange of the two operator sets. We have found it convenient to undo this shift by redefining f − (r), the part of the f-matrix belonging to the anti-symmetric operators, to e iπr f − (r). The resulting combined Fourier spectrum for f + and e iπr f − has strong peaks at k = 2k F and a smaller peak at k = 0, in agreement with the result from [2].
Based on the Fourier spectrum, we rewrite the fitting form (4.10) as with real numbers A µµ > 0 and B µµ , where we expect γ > γ, due to the relative sharpness of the peaks in the Fourier spectrum. The non-linear fitting over the full range of distances is done in several steps to also include the decaying part at long distances on an equal footing. First, the data is rescaled by r +γ , where we obtained γ from a simple power-law fit, in order to be able to fit the oscillations for all distances with comparable accuracy. Then we fit the rescaled data to (6.3), where initially we use the information from the Fourier spectrum in keeping k fixed to k = 2k F , but finally also release the constraint on k. This procedure showed best results, with relative error bounds up to 2%. The uncertainties are largest for the second term in (6.3) as it acts mainly on short distances, having γ > γ.
The results of this fitting procedure are depicted in figure 6, for all 18 nonzero elements of the f-matrix. We see that the leading power-law exponents deviate from the γ (a) 1.25 13 31 54 66 44 55 22 11 45 32 56 46 12 23 64 65 21  0 0.01 0.02 φ µµ is defined such that it is in the interval [−π, π]. The matrix elements have been grouped according to their relative phases φ µµ (separated by the black, dashed line), which clearly indicate cos and sin behaviour for φ µµ = 0 and φ µµ = ± π 2 , respectively. The solid red lines in panels (a) and (b) show the exponent γ 0 and the amplitude A, respectively, from the single fit (6.4).
fit to the r.m.s. net correlations in (3.5) (compare figure 4) by about 5%. The k-vectors from the non-linear fit are close to k = 2k F and deviate by less than 1%. The fit to the sub-leading second term in (6.3) is not reliable, so we do not show the results for γ here, but note that every fit satisfied γ > γ.
Since most of the exponents γ and amplitudes A µµ are of comparable size, we fit the f-matrix elements to a single γ 0 and A (as well as a single γ 0 and B for the second term) for all the f-matrix elements, using the Ansatz: sin(kr) cos(kr) − sin(kr) cos(kr) − sin(kr) cos(kr) sin(kr) cos(kr)   figure 6. Fitting to (6.4) gives an error of about 10%, with largest errors arising for the f-matrix elements where A µµ deviates strongly from A (see figure 6). For the filling ν = 0.286 we find γ 0 = 1.26 and A = 0.06. The values of γ 0 and B are unreliable in that the results from several fittings differ by about 30%, but still it holds that γ 0 > γ 0 .
The form of (6.4) allows us to understand why the r.m.s. net correlations displayed in figure 4 show some residual oscillations, instead of decaying completely smoothly, as anticipated in section 1.2. The reason is that (6.4) contains 10 cos(kr) terms but only 8 sin(kr) terms. Although any two such terms oscillate out of phase, as illustrated in figure 7, the cancellation of oscillations will thus not be complete. Instead, the r.m.s. net correlations contain a factor [8 + 2 cos 2 (kr)] 1 2 (compare to (3.5)), which produces relative oscillations of about 10%, in accord with figure 4. (The fact that the total number of cos(kr) and sin(kr) terms is not equal is to be expected: the total operator Hilbert space per cluster is limited, and its symmetry subspaces might have dimensions not a multiple of 4.) For each pair of wave vectors ±k in each parity sector, the effective operator basis per cluster can be reduced even further, from 3 operators to one conjugate pair of operators. This can be seen by rewriting (6.4) as follows: with the matrices f + and f − defined as Note that both f + and f − are matrices of rank one with eigenvalues 3 2 , 0 and 0. The eigenvectors with eigenvalue 3 2 are 1 √ 3 (1, i, 1) and 1 √ 3 (1, 1, i), respectively. Thus, by transforming to an operator basis in which f ± is diagonal, one finds that in both the even and the odd sector, the dominant correlations are actually carried by only a pair of operators, namely 1 √ 3 (Ô 1 +iÔ 2 +Ô 3 ) and its hermitian conjugate, and 1 √ 3 (Ô 4 +Ô 5 +iÔ 6 ) and its hermitian conjugate, respectively. This result, whose precise form could hardly have been anticipated a priori, is a pleasing illustration of the power of a CDM analysis to uncover nontrivial correlations.

One-particle correlations
The correlations in the 1P sector are exponentially decaying, as already mentioned in section 5.3. The reason for this was given in [1] and is the key to understanding the operators and correlations in this sector. In the limit where the fermions are all paired, the only possible way to annihilate one at x and create one at x > x , such that the initial and final states are both paired, is that every rung in the interval (x, x ) has a fermion (necessarily on alternating legs). These fermions can be grouped as pairs in two different ways: (x, x + 1), (x + 2, x + 3), . . . , (x − 2, x − 1) in the initial state, but (x + 1, x + 2), . . . , (x − 1, x ) in the final state. (Notice this requires that x and x have the same parity.) [1] showed that the probability of such a run of filled sites decays exponentially with its length.
Applying the operator analysis in this sector using the eigenvalue decomposition in (4.6) gives a series of fourfold degenerate eigenvalues for both clusters, see table 2 for cluster A. The table for cluster B is exactly the same. For a specific eigenvalue, also the operators for cluster B (residing at rungs (x , x + 1)) are the same as for cluster A (residing at rungs (x, x + 1)), but with mirrored rungs, i.e. an operator acting on rungs (x, x + 1) acts in the same fashion on rungs (x + 1, x ).
Looking more closely, the first four operators annihilate or create a particle on rungs x + 1 or x , respectively, thereby breaking or regrouping bound pairs residing on (x + 1, x + 2) or (x − 1, x ), respectively. The second set of four operators annihilates or creates a particle on rungs x or x + 1, thereby breaking or regrouping bound pairs residing on rungs (x, x + 1) or (x , x + 1). For a given odd separation x − x, the combination of x + 1 with x requires the smallest number of pairs to be present in between the two clusters. The alternative combination is x with x + 1, which requires an additional pair in between (see figure 8). We could estimate their weights since the relative probability of an extra pair is the factor associated with increasing the separation by two. Since the correlations decay roughly as ∼ 10 −r (see figure 10), we predict two orders of magnitude. Similarly, when x − x is even, we get at mixture of the first and second four operators (see figure 8). This explains the difference in the weights of the two operator sets. Thus, it turns out that for the 1P correlations a cluster size of one rung would already have been large enough to reveal the dominant correlations. We will hence use as operator basiŝ together with their hermitian conjugates. (The fact that our operator basis consists only of operators acting on a single rung implies that it would have been sufficient to use single-rung clusters. However, for the sake of consistency with the rest of our analysis, we retain two-rung clusters here, too.) The f-matrix based on these four operators (per cluster) is diagonal with equal entries for a given distance r. Its Fourier transform (see figure 9) gives a result distinct from the Fourier transform for CD and 2P correlations. The dominant wave vectors are k = ±k F and k = π ± k F , where the latter is the product of an oscillation with k = π and an oscillation with k = ±k F . In total we have an oscillation in the correlations of the form (1 + (−1) r )e ±ik F r , i.e. an oscillation with k = ±k F , and every second term being close to zero. The dominant wave vector k = ±k F i s consistent with the usual behaviour of 1P Green's functions.
The reason for every second term being essentially zero is that the dominant hopping in the system, the correlated hopping, always changes the position of a particle by two rungs, so every second position is omitted. The small but finite value for hopping onto intermediate rungs is related to the finite t /t c = 10 −2 that we use. It results in a second oscillation at k = ±k F located at intermediate rungs, whose relative strength compared to the dominant one is about 10 −2 , which is consistent with the ratio t /t c that we used (see figure 10).
We fit the one independent f-matrix element f µ,µ to an exponential decay of the form Ae −r/r 1 (see figure 10), but apart from this we were not able to fit the exact functional dependence on r, especially the oscillations with k = ±k F . The reason for this is the existence of two oscillations where one is zero on every second rung, and that the data range for which reasonable 1P correlations are still present is too small and thus makes it susceptible to numerical noise. This can be seen already in the Fourier spectrum, where we find relatively broad peaks, as a result of the influence of the exponential envelope and the relatively short distance range available.

Two-particle correlations
The operator subspace for 2P (∆N = 2), in a cluster including two rungs has the comparatively small dimension of four due to the infinite nearest-neighbour repulsion (see figure 2). These areĉ ↑,xĉ↓,x+1 ,ĉ ↓,xĉ↑,x+1 and their hermitian conjugates. In the present case of dominating t c , these operators represent the creation-and annihilationoperators of bound pairs [2]. The operator analysis yields exactly the same four operators with degenerate weight for all distance regimes for both cluster A and B. The four operators are 1/ √ 2 (ĉ ↑,xĉ↓,x+1 ±ĉ ↓,xĉ↑,x+1 ) together with their hermitian conjugates, and they already represent the symmetric and antisymmetric combinations of the operators mentioned above.
The f-matrix (4.7) is diagonal in the basis of the four operators, with equal strength of correlations for a fixed distance apart from a possible sign. This may be expected, given the similar structure of the operators.
As for the CD correlations (∆N = 0), we apply a Fourier transform on the f-matrix (see figure 11) to identify the dominant wave vectors. Again, we find two spectra of similar form but shifted by π with respect to each other. Consequently we redefine f + to e iπr f + , the part of the f-matrix belonging to the symmetric operators. Thus, we obtain one leading peak at k = 0 and sub-leading peaks at k = 2k F . Given the similar structure of the Fourier spectrum to that of the CD correlations, we fit the elements of the f-matrix to the form (6.3), but now expect γ < γ from the relative sharpness of the peaks. Already at the level of the f-matrix elements we find an overall leading decay with residual oscillations, whose relative magnitude becomes smaller at large distances (since γ < γ). Since all matrix elements are the same after redefining f + , it is sufficient to fit |f µ,µ | for a given µ, which will have dominant k-vectors k = 0 and k = ±2k F . The fit has errors of less than 5% throughout, with results as shown in figure 12. The overall behaviour is very similar to the one already found from the r.m.s. net correlations of   this sector (see figure 4), up to the oscillatory part from the second term in (6.3). We see that the oscillations clearly decay more strongly than the actual strength |f µ,µ |, in accord with γ < γ.
In contrast to the CD correlations (see figure 6.1.2), for the 2P correlations we do not find correlations which oscillate with phases shifted by ∆φ = ±π/2 . This may come from the fact that clusters with the size of two rungs have the minimal possible size to capture 2P correlations. The corresponding operator space has dimension four and the four possible operators are very similar in structure. We expect that for larger clusters and hence a larger operator space, we would find correlations which also oscillate out of phase such that their oscillations cancel in the r.m.s. net correlations , in accord with (1.4).

Comparison to previous results
We are now ready to compare our CDM-based results with those obtained in [2] by Cheong and Henley (CH) from fitting simple correlation functions. The latter were computed exactly in [2] for accessible separations after mapping the large t c model onto a hard-core bosonic system, but the functional forms of the r dependencies were inferred from a purely numerical fitting procedure. Overall, our results for the Hamiltonian (2.1) in the strongly correlated hopping regime agree with [2], in that (i) 2P correlations and CD correlations show power-law behaviour, (ii) the 2P correlations dominate at large distances for the fillings we were investigating, (iii) 1P correlations are exponentially decaying and are negligible over all but very short distances, and (iv) the dominating k-vector, for either 2P or CD sectors, is 2k F . However, the power-law exponents obtained from fitting f-matrix elements to (4.10) and summarized in table 3, clearly deviate from the results in [2] by CH. For the CD correlations, in [2] the dependence of γ 0 on the filling ν was given by the exponent γ CH 0 = 1 2 + 5 2 1 2 − ν , from which our results deviate (see figure 4 a,b) by about 25%. Nevertheless, our results for γ 0 agree qualitatively with this prediction, in that we also find γ 0 to decrease linearly with increasing filling.
The 2P correlations deviate more strongly. For the dominant 2P correlations, CH predicted a constant power-law exponent of γ CH 2 = 1 2 independent of filling, coming from a universal correlation exponent for a chain of tightly-bound spinless fermion pairs [8]. In contrast, we obtain a larger exponent (see figure 4 a,b) for given fillings. Our result for γ 2 linearly decreases as the filling gets smaller and appears to approach 1 2 only in the limit ν → 0. We also explicitly calculated the same correlation function as investigated in [2] but found a stronger decay than the r − 1 2 suggested there. We do not know whether the deviation is an artifact of the boundaries of our finite system, or whether the mapping used in [2] to a set of hardcore bosons might have omitted an important contribution.
Moreover, it may be noted that by extrapolating the exponents in a linear fashion towards large fillings (ν → 1 2 ), it appears that for fillings larger than ∼ 0.35 eventually the CD correlations dominate over 2P correlations (see figure 13). This conclusion has also been found in [9] which similarly addresses diatomic real space pairing in the context of superconductivity. Their discussion, however, is not specifically constrained  Figure 13. The power-law exponents for CD correlations (γ 0 , blue symbols) and 2P correlations (γ 2 , red symbols) obtained from the r.m.s. net correlations for several fillings ν. We used chain lengths of N = 100 (circles), N = 150 (crosses), and N = 200 (triangles). The dashed blue and red lines are linear fits to our numerical data for γ 0 and γ 2 , respectively. The solid blue and red lines show the corresponding predictions of Cheong and Henley [2]. For the 2P correlations, our data implies a linear ν-dependence going from 1 2 for ν = 0 to 3 2 for ν = 1 2 . This crossover from 1 2 to 3 2 is predicted by Cheong and Henley as a sub-leading contribution, without giving an explicit functional dependence on ν. The two linear ν-dependencies imply that for large fillings CD correlations should become dominant over 2P correlations. Unfortunately, we do not have been able to obtain reliable data in that regime, because the r.m.s. net correlations showed strong oscillations here, contrary to our expectations from section 1.2.
to one-dimensional systems, and one may wonder how the specific choice of parameters compare. As the filling approaches 0.5 in an excluded-fermion chain, it is appropriate to think about the degrees of freedom as impurity states or holes in the crystalline matrix of pairs [9]. Then the natural length scale is the spacing between holes. The longer that spacing gets (it diverges as ν → 0.5), the larger also the system under investigation must be in order to reach the asymptotic limit. In other words, to see proper scaling behavior in a uniform way, the system size should increase proportional to 1/(0.5 − ν). In our case the data became unreliable for ν 0.4 (see figure 13). On the other hand, for certain fillings ν 0.4, we calculated the power-law exponents for CD and 2P correlations for ladders of length N = 150 and N = 200 (this data is also included in figure 13) and did not find different behaviour compared to out original data for ladders of length N = 100.

Conclusions
Summarizing, we found that the CDM is a useful tool to detect dominant correlations in a quantum lattice system. Starting from a ground state calculated with DMRG, we extracted all the important correlations present in our model system. We developed a method which, first, determines the distance-independent operators on each cluster that carry the dominant correlations of the system, and second, encodes the distancedependence of the correlations in the f-matrix. The latter is then analyzed in terms of decaying and oscillatory terms to extract the long-range behaviour of the correlations.
We saw that the size of the clusters A and B is a limitation of the method as it constrains the analysis to local operators. For some kind of correlations, however, larger clusters are needed to capture the relevant physics. This is not too easily implemented as it requires significantly more resources. As a possible alternative and as an outlook for possible future work, one may think of using a different cluster structure: one cluster as before and one "super-cluster" representing a larger continuous part of the system including one boundary. As MPS introduces, for each site, effective left and right Hilbert spaces describing the part of the chain to the left and to the right of that site, the description of such a super-cluster should be straightforward. The resulting effective density matrix describing a large part of the system can be calculated accordingly.
Overall, DMRG is a suitable method to calculate the CDM. The latter is easily and efficiently calculated within the framework of the MPS. The explicit breaking of (i) translational invariance by using finite system DMRG and (ii) a discrete symmetry of the model, lead us to develop certain strategies to restore these broken symmetries. The smoothing of the boundaries can still be further optimized, or be replaced by periodic boundary conditions. However, we do not expect that this will have significant influence on the conclusions drawn.

A. The variational matrix product state approach
This appendix offers a tutorial introduction to the variational formulation of DMRG for finding the ground state of a one-dimensional quantum lattice model, , based on matrix product states (MPS). It also explains how this approach can be used to efficiently calculate the CDM. We point out all the important properties of the MPS and explain how to perform basic quantum calculations such as evaluating scalar products and expectation values, as well as determining the action of local operators on the MPS and constructing a reduced density matrix. We explain how a given MPS can be optimized in an iterative fashion to find an excellent approximation for the global ground state. We also indicate briefly how the efficiency of the method can be enhanced by using Abelian symmetries.
We would like to emphasize that we make no attempt below at a historical overview of the DMRG approach, or at a complete set of references, since numerous detailed expositions of this approach already exist in the literature (see the excellent review by U. Schollwöck [4]). Our aim is much more modest, namely to describe the strategy implemented in our code in enough detail to be understandable for interested nonexperts.

A.1. Introduction
Quantum many-body systems deal with very large Hilbert spaces even for relatively small system sizes. For example, a one-dimensional quantum chain of N spin 1 2 particles forms a Hilbert space of dimension 2 N , which is exponential in system size. For quantum lattice models in 1D a very efficient numerical method is the density matrix renormalization group (DMRG), introduced by Steven R. White [3]. The problem of large Hilbert space dimension is avoided by an efficient description of the ground state, which discards those parts of the Hilbert space which have negligible weight in the ground state. In this manner the state space dimension of the effective description becomes tractable, and it has been shown that this produces excellent results in many quasi one-dimensional systems.
The algebraic structure of the ground state for one-dimensional systems calculated with DMRG is described in terms of matrix product states (MPS) [10,11,12,5,13]. The origin of this MPS structure can be understood as follows (a detailed description will follow later): pick any specific site of the quantum lattice model, say site k, representing a local degree of freedom whose possible values are labeled by an index σ k (e.g., for a chain of spinless fermions, σ k = 0 or 1 would represent an empty or occupied site). Any many-body state |ψ of the full chain can be expressed in the form where |l k and |r k are sets of states (say N l and N r in number) describing the parts of the chain to the left and right of current site k, respectively, and for each σ k , l k r k and dimension N l × N r . Since such a description is possible for any site k, the state |ψ can be specified in terms of the set of all matrices A [σ k ] , resulting in a matrix product state of the form One may now seek to minimize the ground state energy within the space of all MPS, treating the matrix elements of the A-matrices as variational parameters to minimize the expectation value ψ| H |ψ . If this is done by sequentially stepping through all matrices in the MPS and optimizing one matrix at a time (while keeping the other matrices fixed), the resulting procedure is equivalent to a strictly variational minimization of the ground state energy within the space of all MPS of the form (A.2) [5,10,11,12,13]. If instead the optimization is performed for two adjacent matrices at a time, the resulting (quasi-variational) procedure is equivalent to White's original formulation of DMRG [5,10,11,12,13]. The MPS based formulation of this strategy has proven to be very enlightening and fruitful, in particular also in conjunction with concepts from quantum information theory [5].
In general, such an approach works for both bosonic and fermionic systems. However, to be efficient the method needs a local Hilbert space with finite and small dimension, limiting its applicability to cases where the local Hilbert space is finite dimensional a priori (e.g. fermions or hard-core bosons) or effectively reduced to a finite dimension, e.g. by interactions. For example, such a reduction is possible if there is a large repulsion between bosons on the same site such that only a few states with small occupation number will actually take part in the ground state. For fermions, on the other hand, the fermionic sign must be properly taken care of. The anti-commutation rules of fermionic creation and annihilation operators causes the action of an operator on a single site to be non-local because the occupations of the other sites have to be accounted for. To simplify the problem, a Jordan-Wigner transformation [14] can be used to transform fermionic creation and annihilation operators to new operators that obey bosonic commutation relations for any two operators referring to different sites. This greatly simplifies the numerical treatment of these operators as fermionic signs can be (almost) ignored.
Before outlining in more detail the above-mentioned optimization scheme for determining the ground state (see section A.3), we present in section A.2 various technical ingredients needed when working with MPS. i.e. the matrix decomposition is a 1 × 1 matrix which is a scalar. If these A-matrices are sufficiently large this decomposition is formally exact, but since that would require A-matrices of exponentially large size, such an exact description is of academic interest only. The reason why the A-matrices are introduced is that they offer a very intuitive strategy for reducing the numerical resources needed to describe a given quantum state. This strategy involves limiting the dimensions of these matrices by systematically using singular-value decomposition and retaining only the set of largest singular values. The A-matrices can be chosen much smaller while still giving a very good approximation of the state |ψ .
Selecting a certain site k, the state can be rewritten in the form (A.1). The  This reduces the resources used to describe a state from O(d N ) for the full many-body Hilbert space down to O(N D 2 d). This is linear in the system size, assuming that the size required for D to accurately describe the state grows significantly slower than linearly in N . This, in fact, turns out to be the case for ground state calculations [15]. Details of this truncation procedure and estimates of the resulting error are described in section A.2.5.

A.2.2. Global view and local view
Matrix product states can be viewed in two alternative ways: a global view and a local view. Both views are equivalent and both have their applications. In the global view the state is expressed as in (A.2), i.e. the effective Hilbert spaces have been used 'only' to reduce resources. The state is stored in the A-matrices, but the effective basis sets will be contracted out. This perception has to be handled very careful, because contracting out the effective basis sets leads to higher costs in resources! In the local view the state is expressed as in (A.1). It is called local because there is one special site, the current site, and all other sites are combined in effective orthonormalized basis sets. Usually, the local view is used iteratively for every site. In this perception, we need effective descriptions of operators contributing to the Hamiltonian acting on other sites than the current site (see section A.2.8).

A.2.3. Details of the A-matrices
The A-matrices have some useful properties that hold independently of the truncation scheme used to limit the effective Hilbert spaces. First of all, we notice that by construction dim(H r k−1 ) ≡ dim(H l k ), otherwise the matrix products in (A.2) would be ill defined. Based on this, we can find another interpretation of the A-matrices in the local view. The part of the chain to the left of site k (where k is far from the ends for simplicity) is described by the effective basis |l k , which is built of truncated A-matrices: The A [σ k−1 ] -matrix maps the effective left basis |l k−1 together with the local |σ k−1 basis onto the effective left basis |l k ! The same argument applied on the effective right basis of site k leads to the transformation of |r k+1 and |σ k+1 onto |r k via the A [σ k+1 ] -matrix: So far, this may be any transformation, but in order to deal with properly orthonormal basis sets, we may impose unitarity on the transformation (see below). The A-matrices towards the ends of the chain have to be discussed separately. The use of open boundary conditions implies that we have a 1-dimensional effective state space to the left of site one and the right of site N , respectively, both representing the empty state. This implies that dim(H l 1 ) = 1 = dim(H r N ). Moving inwards from the ends of the chain, the effective Hilbert spaces acquire dimension d 1 , d 2 , . . . until they become larger than D and need to be truncated. Correspondingly, the dimension of matrix . In these cases we simply choose A (l k σ k )r k = 1 and A l k (r k σ k ) = 1, respectively.
Summarizing, the A-matrices have two functions. If site i is the current site in (A.1), the A [σ i ] -matrices represent the state, i.e. its coefficients specify the linear combination of basis states |l k , |σ k and |r k . On the other hand, if not the current site, the Amatrices are used as a mapping to build the effective orthonormal basis for the current site, as we describe next: Orthonormal basis sets In the local view, the whole system is described by the Amatrices of the current site k in the effective left basis, the effective right basis, and the local basis of site k. A priori, the basis states form an orthonormal set only for the local basis set, but we may ask for the effective basis sets |l and |r ‡ to be orthonormal, too, i.e. require them to obey: This immediately implies the following condition on the A [σ j ] -matrices, using (A.4) and (A.5) (for a derivation, see section A.5.1): The orthonormality (A.6) for both the left-and right basis states holds only for the current site. For the other sites there is always only one orthonormal effective basis.
Graphical representation Matrix product states can be depicted in a convenient graphical representation (see figure A2). In this representation, A-matrices are displayed as boxes and A [σ k ] is replaced by A k for brevity. Indices correspond to links from the boxes. The left link connects to the effective left basis, the right link to the right one, and the link at the bottom to the local basis. Sometimes indices are explicitly written on the links to emphasize the structure of the sketch. Connected links denote a summation over the indices (also called contraction) of the corresponding A [σ] -matrices. At the boundaries of the chain, a cross is used to indicate the vacuum state. A.2.4. Orthonormalization of effective basis states We now describe how an arbitrary MPS state can be rewritten into a form where its local view with respect to a given site has orthonormal left-and right basis states. It should be emphasized that this really just amounts to a reshuffling of information among the state's A-matrices without changing the state itself, by exploiting the freedom that we always can insert any X −1 X = 1 at any position in the matrix product state without altering it. ‡ From now on the index k is only displayed when several sites are involved. For the current site or in the case when only one A-matrix is considered the index will be dropped.
Assume site k to be the current site and assume that it has an orthonormal left basis (the latter is automatically fulfilled for k = 1). We need a procedure to ensure that, when the current site is switched to site k+1, this site, too, will have an orthonormal left basis. (This is required for the orthonormality properties used in the proof in section A.5.1. A similar procedure can be used to ensure that site k − 1 has an orthonormal right basis provided k has such a basis.) For this purpose we use the singular value decomposition (SVD, see section A.5.2) for which we have to rewrite A [σ k ] l k r k by fusing the indices l k and σ k : where m, n and r k have the same index range (see figure A3). Specifically, u fulfills which is equivalent to the orthonormality condition (A.7) for the A [σ k ] -matrices. Figure A3. Singular value decomposition of the A-matrices

SVD
As u replaces A [σ k ] and sv † is contracted onto A [σ k+1 ] , this leaves the overall state unchanged (for a graphical depiction see figure A4): Figure A4. Rearrangement of the A-matrices to switch the current site from site k to k + 1.
Site k +1 now has an orthonormal effective left basis. A similar procedure works for the effective right basis, see figure A5. To obtain an orthonormal effective left basis for the current site k, we start with the first site, update A [σ 1 ] and A [σ 2 ] , move to the next site, update A [σ 2 ] and A [σ 3 ] , and so on until site k − 1. For an orthonormal effective right basis, we start from site N and apply an analogous procedure in the other direction.
If the state |ψ is in the local description of site k with orthonormal basis sets |l k , |σ k and |r k , it is now very easy to change the current site to site k ± 1, with corresponding new orthonormal basis sets |l k±1 , |σ k±1 , |r k±1 . Suppose we want to change the current site from site k to site k + 1. Following the procedure described above, site k + 1 already has an orthonormal right basis and all sites left of site k fulfill the orthonormality condition. All that is left to do, is to update site k and k + 1 to obtain an orthonormal left basis for site k + 1. This is called a switch of the current site from site k to k + 1. The switch from site k to site k − 1 is done analogously.
A.2.5. Hilbert space truncation A central ingredient in the variational optimization of the ground state (see section A.3.1 below) is the truncation of the effective Hilbert spaces associated with a given A-matrix. The strategy for truncating the effective Hilbert spaces is completely analogous to the original DMRG formulation [11]. The DMRG truncation scheme is based on discarding that part of the Hilbert space on which a certain density matrix has sufficiently small weight. There are two ways to obtain an appropriate reduced density matrix: two-site DMRG [3,4] and one-site DMRG [4]. The crucial difference between the two is that one-site DMRG is strictly variational in the sense that the energy is monotonically decreasing with each step,, whereas in two-site DMRG the energy may (slightly) increase in some steps, but with the advantage that the cutoff dimension can be chosen dynamically in each step.
Two-site DMRG Two-site DMRG arises when variationally optimizing two sites at a time. We consider two current sites, say k and k + 1, and we may choose the cutoff dimension site-dependent: D → D k ≡ dim(H l k ). Following section A.2.4, we assume site k to have an orthonormal left basis and site k + 1 to have an orthonormal right basis. After contracting the indices connecting A [σ k ] and A [σ k+1 ] (see figure A6), the state is described by A [σ k σ k+1 ] l k r k+1 . In this description we may optimize the ground state locally by variationally minimizing the ground state energy with respect to A again. This can be accomplished via singular value decomposition (see section A.5.2) by fusing the indices l k , σ k → (l k σ k ) and r k+1 , σ k+1 → (r k+1 σ k+1 ) (see figure A6) to obtain A , where i = 1 . . . min(dD k , dD k+2 ). Using the column unitarity of u and the row unitarity of v † (see section A.5.2), we rewrite the state as where the new set of basis states |l i and |r i is orthonormal with l i |l i = δ i i and r i |r i = δ i i . This representation of the state may be seen as residing on the bond between k and k + 1, with effective orthonormal basis sets for the parts of the system to the left and right of the bond. Reduced density matrices for these parts of the system, obtained by tracing out the respective complementary part, have the form: The standard DMRG truncation scheme amounts to truncating ρ [L] and ρ [R] according to their singular values s i . We could either keep all singular values greater than a certain cutoff, thereby specifying a value for D k+1 between 1 and min (dD k , dD k+2 ), or alternatively choose D k = D to be site-independent for simplicity. This step makes the method not strictly variational, since we discard some part of the Hilbert space which could increase the energy. It turns out that this potential increase of energy is negligible in practice. We can obtain a measure for the information lost due to truncation by using the von Neumann entropy S = − tr (ρ ln ρ), given by where s 2 i = 1 due to the normalization of |ψ . Figure A6. Procedure for site update within two-site DMRG. The grey line under the s indicates that s is the diagonal matrix of singular values.

SVD
One-site DMRG One-site DMRG arises when variationally optimizing one site at a time. In contrast to two-site DMRG, one-site DMRG does not easily allow for dynamical truncation during the calculation. (It is possible in principle to implement the latter, but if one decides to use dynamical truncation, it would be advisable to do so using twosite DMRG.) The truncation is fixed by the initial choice of D, but it is still possible to determine an estimate on the error of this truncation by analyzing the reduced density matrix. Starting from an expression for the full density matrix in the local view (current site k with orthonormal effective basis sets) lr |l |σ |r lr A [σ ] l r * |l l | |σ σ | |r r | , (A.14) we trace out the effective right basis and obtain a reduced density matrix for the current site and the left part of the system: This reduced density matrix carries the label l k+1 because it corresponds precisely to the density matrix |l k+1 l k+1 . So if we switch the current site from site k to site k + 1, we can check the error of the truncation of H l k+1 . Fusing the indices l and σ, we obtain We do not need to diagonalize the coefficient matrix AA † to obtain the largest weights in the density matrix, because we get its eigenvalues as a byproduct of the following manipulations anyway [4]. To switch the current site we need to apply a singular value decomposition (see section A.2.4) and obtain A = usv † (this is not the usual Amatrix, but the index-fused form). This directly yields AA † = usv † vsu † = us 2 u † , which corresponds to the diagonalization of ρ [l k+1 ] , implying that the weights of the density matrix are equal to s 2 . Of course this works also for the right effective basis. With such an expression, we can check whether the effective Hilbert space dimension D of H l k+1 is too small or not. For example, we could ask for the smallest singular value s D to be at least n orders of magnitude smaller than the largest one s 1 , i.e. the respective weights in the density matrix would be 2n orders of magnitude apart. If the singular values do not decrease that rapidly, we have to choose a greater D.
A.2.6. Scalar product The scalar product of two states |ψ and |ψ is one of the simplest operations we can perform with matrix product states. It is calculated most conveniently in the global view because then we do not need to care about orthonormalization of the A-matrices: using the orthonormality of the local basis σ k |σ l = δ kl δ σ k σ k . In principle the order in which these contractions are carried out is irrelevant, but in practice it is possible to choose an order in which this summation over the full Hilbert space is carried out very efficiently by exploiting the one-dimensional structure of the matrix product state (see figure A7 for a graphical explanation). For details on the numerical costs, see section A.5.3. In method (a), after contracting all A-matrices of |ψ and |ψ , we have to perform a contraction over the full Hilbert space, i.e. a 1 × d N matrix is multiplied with a d N ×1 matrix. This contraction is of order O d N , which is completely unfeasible for practical purposes. In method (b) the most 'expensive' contraction is in the middle of the chain, say at site k, and it is of order O (dD 3 ). Here the A-matrices are viewed as three-index objects A l k r k σ k with dimension D × D × d. All sites left of site k are represented by a D × D matrix, say L l k l k . Contracting this with the matrix at site k yields the object l k L l k l k A l k r k σ k , which has dimensions D × D × d, and since the sum contains D terms, the overall cost is O (dD 3 ). Thus, in practice, method (b) is rather efficient and renders such calculations feasible in practice. Partial product Sometimes it is required to calculate a product over only a part of the matrix product state. This is done the same way as the scalar product Notice that P [L k ] and P [R k ] are matrices in the indices l k and r k , respectively (see figure A8). In fact, they correspond to the overlap matrices l k |l k and r k |r k , respectively. Figure A8. Partial products associated with site k.
A.2.7. Reduced density matrix The pure density matrix given by the matrix product state |ψ is defined as ρ = |ψ ψ|. To describe only a part of the system, we need to calculate the reduced density matrix. Let I be a set of sites and σ s = {σ k∈I } a fused index for their local states. Tracing out all other sites with combined index σ b = {σ k / ∈I } we obtain This is a completely general expression, but in the cases where I = {k} or I = {k, k } it reduces to (see figure A9) A similar strategy can be used to calculate the density matrices needed for the main text, by contracting out the σ k 's for all sites except those involved in the clusters A, B or A ∪ B. In fact, (A.23) givesρ A∪B for two clusters of size one at sites k and k .
where the only condition to derive these results, was that site k − 1 has an orthonormal effective left basis. Similarly, if the (k −1)-left-representation of an operator C is known, its k-left-representation can be obtained via (see figure A10) Equation (A.24) and (A.25) can be used iteratively to transcribe the i-localrepresentation of B into its k-left-representation for any k > i (see figure A11). This  reasoning also applies to the right site of site k and so it is possible to obtain a description of any local operator on any site. To obtain a description of a pair of local operators acting on different sites, we have to transcribe them step by step. Let site k be the current site with orthonormal effective basis sets and B, C two operators acting locally on site i and j respectively (i < j < k). First we obtain the j-left-representation of B, namely B l j l j , as described above. Then both operators are transformed together into the (j + 1)-left-representation (see figure A12), which in turn can be transformed iteratively into the desired k-left-representation of the operators B and C. Figure A12. The (j + 1)-left-representation of the operators C, given in the j-localrepresentation, and B, given in the j-left-representation.
A.2.9. Local operators acting on |ψ Any combination of operators can be calculated directly in the global view or in the local view via the effective descriptions introduced in the previous section.

Global view
The operators, known in the local basis of the site they are acting on, are contracted directly with the corresponding A-matrix. For example, the formula for a nearest neighbour hopping term c † k c k+1 (see figure A13) reads as Local view Let k be the current site with orthonormal effective basis sets. If we want to evaluate operators acting on other sites than the current site k, we need an effective description of these operators in one of the effective basis sets of site k to contract these operators with the A-matrix of the current site. For example, to calculate the action of the nearest neighbour hopping term c † k c k+1 on |ψ = A lr |l |σ k |r , we need (c † k ) σ k σ k and (c k+1 ) r r to obtain (see figure A13) Figure A13. The nearest neighbour hopping term c † k c k+1 acting on |ψ in (a) the global view and (b) the local view.
A.2.10. Expectation values Expectation values are merely the scalar product between the state with itself including the action of an operator and can be easily worked out in both the global and the local view (see figure A14). Since both methods are equivalent, the local variant is much more efficient as it involves much less matrix multiplications. However, it requires careful orthonormalization of the remainder of the A-matrices. The iterative scheme, introduced in section A.3, allows for that and works in the local picture.

A.3. Variational optimization scheme
The basic techniques introduced in the previous sections are the building blocks for DMRG sweeps, an iterative scheme to determine the ground state in the usual DMRG sense. This scheme starts at some site as current site, for example the first site where truncation occurs, and minimizes the energy of |ψ with respect to that site. Afterwards the current site is shifted to the next site, and the energy of |ψ with respect to that site is minimized. This is repeated until the last site where truncation occurs is reached and the direction of the switches is reversed. When the starting site is reached again, one sweep has been finished (see figure A15). These sweeps are repeated until |ψ converges. A.3.1. Energy minimization of the current site In order to find the ground state of the system we have to minimize the energy E = ψ| H |ψ of the matrix product state |ψ with the constraint that the norm of |ψ must not change. Introducing λ as Lagrange multiplier to ensure proper normalization, we arrive at the problem of determining min |ψ ( ψ| H |ψ − λ ψ|ψ ) . (A.29) In the sweeping procedure introduced above, the current site is changed from one site to the next and the energy is minimized in each local description. where H l r σ lrσ = l | σ | r | H |l |σ |r is the Hamiltonian expressed in the two orthonormal effective basis sets and the local basis of the current site. The multidimensional minimization problem (A.29) has been transformed to a local minimization problem where one A-matrix (or two) is optimized at a time and all others are kept constant. Such a procedure could, in principle, cause the system to get stuck in a local minimum in energy, but experience shows that the procedure works well [4], especially in the presence of a gap. The matrix elements H l r σ lrσ may be calculated easily using the techniques introduced in section A.2 (see section A.3.2 for details). Changing to matrix notation and replacing λ with E 0 in anticipation of its interpretation as an energy, we obtain an eigenvalue equation: lr |l |σ |r . (A.32) The minimization problem reduces to a local eigenvalue problem, which can be solved by standard techniques. The full Hilbert space of the current site has dimension dD 2 and may become large, but it is not necessary to determine the full spectrum of H, since we are interested only in the ground state. The Lanczos algorithm is an effective algorithm to achieve exactly that. The advantage of this algorithm is that we only have to compute H |ψ , which saves much effort. The Lanczos algorithm produces as output the ground state eigenvalue and eigenvector. The latter gives the desired optimized version of the matrix A σ lr , which then has to be rewritten (with or without Hilbert space truncation, as needed) into a form that satisfies the orthonormality requirements of the left and right basis sets, as described in section A.2.4.
A.3.2. Sweeping details Before the actual sweeping may be started we have to set up an initial state, prepare a current site with orthonormal effective basis sets and calculate effective descriptions of operators which are part of the Hamiltonian. After this initialization we may determine the ground state with respect to this current site and shift the current site to the next site. That current site again has orthonormal effective basis sets due to the switching procedure introduced in section A.2.4, but we also need effective representations of the operators acting in the Hamiltonian. At this step the structure of the matrix product state saves much effort, as most of the needed representations are already calculated.
Structure of the Hamiltonian terms The Hamiltonian H l r σ lrσ , acting in the space spanned by the states |l , |σ , |r , breaks up into several terms: where the indices denote on which parts of the system the respective term acts on (L and R indicate left and right of the current site, respectively, • indicates action on the current site). Of course, the six terms of (A.33) depend on the current site k: H The terms (H L ) l l and (H R ) r r contain all terms which involve only sites k < k and k > k, respectively. The iterative structure of the method directly yields the following equalities: where the terms on the right hand side are meant to be expressed in the effective basis of the operator on the left hand site (see figure A17). Initialization First of all we need an initial matrix product state, which is most conveniently chosen to consist of identity transformations at the ends of the chain (see section A.2.3) and random A-matrices for the rest of the chain. We take the first site where Hilbert space truncation is applied as current site k and obtain an orthonormal effective right basis (the effective left basis is already orthonormal) using the orthonormalization procedure introduced in section A.2.4 starting from site N . Additionally it is convenient, while dealing with site N , to calculate and store the operator H This ensures, when the sweeping procedure reaches site N − 1, that all necessary operators are already calculated. This is repeated from site N down to site k + 1, and similarly for the sites k < k in the other direction. The result of these initialization steps is that we have a current site k with orthonormal effective basis sets, effective descriptions of the Hamiltonian terms H  Complete ground state calculation The methods introduced above make the procedure to determine the ground state very efficient as the global problem is mapped onto many local problems involving only a few terms to calculate. The iterative structure of the matrix product states and the effective Hamiltonian terms strongly increase the efficiency. A full ground state calculation consists of: (i) Initialization as described above (ii) Full sweeps from site K to site K and back to site K, with sites K and K the first and last site where the effective Hilbert spaces are truncated.
(iii) After each sweep i the overlap ψ i−1 |ψ i between the state before and after the sweep is calculated. If the matrix product state does not change any more, stop the sweeping. A criterion, for example, for when to stop would be to require that where is a small control parameter, typically of order 10 −10 .
Numerical costs The step with the most impact on the numerical costs of the algorithm is the calculation of H |ψ in the Lanczos method. This method is an iterative scheme using several Lanczos steps, of which usually less than 100 are needed for one ground state calculation. Each Lanczos step calculates H |ψ exactly once. This calculation basically consists of elementary matrix multiplications, see section A.5.3 for details on the numerical costs of such calculations. The six terms introduced in (A.33) are not all equally time consuming. Most of them contain identity maps which do not need to be carried out and thus the term H L•R is the most time consuming, requiring operations of order O(dD 2 (2D + d)). The total numerical cost for the minimization process is where N Sweep is the number of sweeps, N the chain length and N Lanczos the number of Lanczos steps. In practice the cutoff dimension is significantly higher than the local Hilbert space dimension d and thus (A.37) is nearly linear in d.

A.4. Abelian symmetries
Matrix product states can be easily adapted to properly account for conserved quantum numbers, representing the global symmetries of the Hamiltonian. We will limit ourselves to Abelian symmetries, meaning that the irreducible representation of the symmetry group is Abelian, as these are easily implemented, which is not necessarily the case for non-Abelian symmetries [16]. An Abelian symmetry allows a quantum number Q to be attached to every state. The property that the symmetry is Abelian manifests itself in that this quantum number is strictly additive. For two states |Q 1 and |Q 2 , the quantum number of the direct product of these two states is given by |Q 1 ⊗ |Q 2 = |Q 1 + Q 2 . For example, if the Hamiltonian commutes with the number operator for the full system, the quantum number Q could represent particle number.
For matrix product states, the introduction of Abelian symmetries has the consequence that the A-matrix A [σ] lr may be written as (A Qσ Q l Qr ) γσ α l βr . Here Q σ , Q l , Q r are the quantum numbers attached to the local, left effective and right effective basis, respectively. The index α l distinguishes different states |Q l , α l characterized by the same quantum number Q l , and similarly for |Q r , β r and |Q σ , γ σ . If A describes, for example, the mapping of the |l -basis of the left block together with the local basis to a combined (truncated) |r -basis, then the only non-zero blocks of the A-matrix are those for which Q σ + Q l = Q r . For the current site, the total symmetry Q tot of the full quantum many-body state manifests itself in that the corresponding A-matrix fulfills Q l + Q r + Q σ = Q tot .
For the handling of matrix product states quantum numbers imply a significant amount of bookkeeping, i.e. for every coefficient block we have to store its quantum number. The benefit is that we can deal with large effective state spaces at reasonable numerical cost. The Lanczos algorithm, in particular, takes advantage of the block structure.
Of course, the treatment of Abelian symmetries is generic and not limited to only one symmetry. We may incorporate as many symmetries as exist for a given Hamiltonian, by writing Q as a vector of the corresponding quantum numbers.
A.5. Additional details A.5.1. Derivation of the orthonormality condition The orthonormality condition (A.7) is easily derived by induction. The starting point is condition (A.6) and we limit to the derivation for the left basis. The derivation for the right basis is analogous.
The induction argument can be initialized with site k = 1 because its effective left basis is already orthonormal as it consists only of the vacuum state. Now, consider the case that site k has an orthonormal effective left basis and construct the condition for site k + 1 to have an orthonormal effective left basis: Condition (A.7) follows with l k+1 l k+1 ! = δ l k+1 l k+1 .
A.5.2. Singular value decomposition The singular value decomposition can be seen as a generalization of the spectral theorem, i.e. of the eigenvalue decomposition. It is valid for any real or complex m × n rectangular matrix. Let M be such a matrix, then it can be written in a singular value decomposition where U is a m × m unitary matrix, S a m × n matrix with real, nonnegative entries on the diagonal and zeros off the diagonal, and V a n × n unitary matrix. The numbers on the diagonal of S are called singular values, and there are p = min (n, m) of them. The singular values are unique, but U and V are not, in general. It is convenient to truncate and reorder these matrices in such a fashion that their dimension are m × p for U , p × p for S (with the singular values ordered in a non-increasing fashion) and n × p for V (i.e. p×n for V † ). A consequence of this truncation is that U or V is no longer quadratic and unitarity is not defined for such matrices. This property is replaced by column unitarity (orthonormal columns) of U and row unitarity (orthonormal rows) for V † -no matter which one is no longer quadratic. In this article all singular value decompositions are understood to be ordered in this fashion.

A.5.3.
Numerical costs of index contractions The numerical costs of matrix multiplications and index contractions of multi-index objects depend on the dimension of both the resulting object and of the contracted indices. In the case of matrix multiplications this is quite simple. Consider a n × m matrix M 1 multiplied by a m × p matrix M 2 . The result is a n × p matrix M : Thus for every entry of M , p 1 times p 2 multiplications have to be done, so that the process is of order O ((p 3 . . . p n ) (p 1 p 2 ) (q 3 . . . q m )).