Bundled matrix product states represent low-energy excitations faithfully
Thomas E. Baker Negar Seif

















theoremTheorem
definition
definitionDefinition

remark
*remarkRemark

corollaryCorollary[theorem]
lemma[theorem]Lemma
proposition[theorem]Proposition

























ε̅

r
N_e
N
d
N_s
N_ML
d 
d' 











Julia
  morekeywords=abstract,break,case,catch,const,continue,do,else,elseif,
      end,export,false,for,function,immutable,import,importall,if,in,
      macro,module,otherwise,quote,return,switch,true,try,type,typealias,
      using,while,
   sensitive=true,
   alsoother=,
   morecomment=[l]#,
   morecomment=[n]#==#,
   morestring=[s]"",
   morestring=[m]'',
[keywords,comments,strings]
    language         = Julia,
    basicstyle       = ,
    keywordstyle     = ,
    numbers=left,
    stringstyle      = ,
    commentstyle     = ,
    showstringspaces = false,
    frame = single, 
FIGs/
[Please direct correspondence to: ]bakerte@uvic.caDepartment of Physics & Astronomy, University of Victoria, Victoria, British Columbia V8P 5C2, CanadaDepartment of Chemistry, University of Victoria, Victoria, British Columbia V8P 5C2, CanadaCentre for Advanced Materials & Related Technologies, University of Victoria, Victoria, British Columbia V8P 5C2, CanadaDepartment of Physics & Astronomy, University of Victoria, Victoria, British Columbia V8P 5C2, Canada§ ABSTRACT 

We consider a set of density matrices. All of which are written in the same orbital basis, but the orbital basis size is less than the total Hilbert space size. We ask how each density matrix is related to each of the others by establishing a norm between density matrices based on the truncation error in a partial trace for a small set of orbitals. We find that states with large energy differences must have large differences in their density matrices. Small energy differences are divided into two groups, one where two density matrices have small differences and another where they are very different, as is the case of symmetry. We extend these ideas to a bundle of matrix product states and show that bond dimension of the wavefunction ansatz for two states with large energy differences are larger. Meanwhile, low energy differences can have nearly the same bond dimensions for similar states.

Bundled matrix product states represent low-energy excitations faithfully
    Negar Seif
    September 9, 2024
=========================================================================



§ INTRODUCTION


Density matrices represent one of the core objects in quantum mechanics. They store a wealth of information about the system and can be useful for solving problems. It is well established that the trace of the density matrix multiplied onto any operator gives the expectation value of the operator for a given state that the density matrix represents <cit.>.

When diagonalized with an eigenvalue decomposition, density matrices are decomposed into a diagonal matrix that contains the orbital occupations of the natural orbitals. The natural orbitals themselves are the eigenvectors of the density matrix. It was originally pointed out by Löwdin <cit.> that natural orbitals were a rapidly converging basis set ( i.e. the lowest eigenvalue converges faster than other choices of a basis with increasing numbers of natural orbitals). Because the solution of natural orbitals requires a ground state wavefunction, it is often computationally costly to obtain them before a computation. So, another basis set is often used.

However, describing the density matrix with a number of states equal to the total Hilbert space size is computationally costly. Reducing the number of degrees of freedom while maintaining accuracy on the result is the main challenge of computational chemistry and solutions of quantum problems on the classical computer in general. This is the foundational idea behind renormalization.

What is considered less often is how natural orbitals for one density matrix describe well or do not describe at all another density matrix. Our goal in this paper is to determine how density matrices when summed over some of the basis functions in a given basis can describe a system. We consider here the idea of bundling together different density matrices. We further consider how accurate those density matrices are if a common basis is used to write both density matrices in full. For example, if a set of m orbtials that have the highest occupation for one density matrix are used to express another density matrix, how accurate can the second density matrix be and what is the best way to characterize it?

The fundamental question that is being asked here is how best can one relate what we will call a bundled set of density matrices, defined as follows.


A bundle of density matrices is defined as a set of density matrices that are all written in the same basis.


This is not the same as an ensemble of states contained in the density matrix.

The fundamental quantity that we want to investigate is whether a notion of closeness ( i.e., a norm) can be defined for the independent density matrices. The result used will be to establish a relationship between the truncation error and the metric distance between two density matrices. This will also be related to the energy difference between two states. The argumentation applies to any local Hamiltonian, which is reasonable for physical systems. We then extend the outcomes of those answers to matrix product states (MPS) to understand how the bond dimension of a bundle of MPSs will behave. This will explain why the bond dimension of the bundled MPS was not explosively large when an algorithm was formulated to solve for excitations in a quantum model <cit.>.

The analysis tools used to formulate the truncation of the bundled MPS apply in principle to bundles of any type of density matrices so long as the eigenstates satisfy the area law of entanglement.  However, we choose to focus on bundles of eigenstates because they apply most readily in entanglement renormalization algorithms. In the following, we use theorems only when they are most relevant to the main thesis statement of the paper. We use definitions throughout to clearly define core concepts and keywords. Corollaries are used to communicate small extensions of the core theorems. Lemmas are used when heavy reliance on results outside of the paper are required to prove the statement and also when those statements are required for proofs later on. 



§ BACKGROUND ON DENSITY MATRICES


The class of problems that we wish to solve are based on the definition of a Hamiltonian that is a self-adjoint operator <cit.>, H, composed of complex coefficients (H∈ℂ^M× M for an M sized Hilbert space). The eigensolutions, ψ, of this operator have the relationship Hψ=Eψ for eigenvalues (energies) E∈ℝ and eigenvectors (wavefunctions) ψ∈ℒ^2 (square integrable) and contain complex numbers <cit.>. The Hamiltonian can be represented by a number of site indices i,j,k,ℓ,…. For many-body Hamiltonians, there is a quartic term (4-indices required) that appears to account for the electron-electron interaction, although the results we derive here will apply for any interaction. We choose to start from many-body Hamiltonians since this will recover a wide class a non-relativistic phenomena that we are interested in.

There is no consideration for divergences ( i.e., points where the evaluation of a quantity is infinite) in any terms as this analysis is solely concerned with models implemented to a finite numerical precision. Thus, all singularities in any interactions are regularized by finite difference approximations.



 §.§ Density matrices


A density matrix can have several connotations. We explicitly define several that are useful here. The type of density matrix that will be used here is the one-body reduced density matrix, although we describe the more general case in many places.


For a given Hamiltonian H, the full density matrix of a system is defined by

    ρ=∑_kη_k|ψ_k⟩⟨ψ_k|

for the kth excitation of the system and some occupation used throughout as 0≤η_k≤1. When projected onto a real space lattice (or other basis) by through a resolution of the identity, the density matrix then assumes the form,

    ρ=   ∑_ ijk𝒜ℬη_k|i𝒜⟩⟨ i𝒜|ψ_k⟩⟨ψ_k|jℬ⟩⟨ jℬ|
    
    =   ∑_ij𝒜ℬρ_ij𝒜ℬ|i𝒜⟩⟨ jℬ|

after defining ρ_ij𝒜ℬ=∑_kη_k⟨ i𝒜|ψ_k⟩⟨ψ_k|jℬ⟩. When 𝒜 and ℬ are a null set, when either variable contains no indices, the density matrix is represented as a one-body reduced matrix ( i.e., only requiring the indices i and j),

    ρ=∑_ijρ_ij|i⟩⟨ j|.

Had more indices been kept, then a higher order density matrix would be represented ( i.e., i,j,k,ℓ for the two-body reduced density matrix). In second quantization, one can simply compute ρ_ij=⟨ c^†_iσ c_jσ⟩ (see Apx. <ref>) in a fermion model with spin σ. The density matrix defined in Def. <ref> is generally representing a  mixed density matrix. 


There are a few types of density matrices that form special types of the above definition. 


In the special case where  ρ^2=ρ or that Tr(ρ^2)=1 the density matrix is called a pure density matrix. Pure states can be represented as ρ=|ψ⟩⟨ψ| for some state ψ.


In quantum chemistry, there is a different normalization convention. The trace of the density matrix is not always one, but instead is the number of particles, N_e, with a particular spin. We find it occasionally useful to refer to an ensemble density matrix that we define as the following. 


Any density matrix whose trace is not one (η_k can assume any positive value).


Note that when the normalization of the ensemble density matrix has a trace that is not 1 that one can have Tr(ρ)≠1 and ρ^2=ρ. For example, three electrons of the same spin in a pure state can be represented in an ensemble density matrix and have Tr(ρ)=3.

In many contexts, the mixed density matrix and the ensemble density matrix are used synonymously. A density matrix in the mixed representation can be thought of as a linear superposition of other states, so the concepts are the same. 

We define explicitly here the ensemble density matrix not only to distinguish between normalization factors in quantum information and quantum chemistry, but we also take this opportunity to highlight that the bundled density matrix from Def. <ref> is not an ensemble density matrix. Specifically, the bundle of density matrices can contain ensemble density matrices, although the density matrices in the bundle can be of any type. 


The diagonalization of a density matrix yields a set of eigenvectors known as the  natural orbitals. For the one-body reduced density matrix, these orbitals represent the density matrix as

    ρ=∑_kε_k|Φ_k⟩⟨Φ_k|

with an eigenvalue ε_k sometimes called an  occupational weight.


At first glance, Eq. (<ref>) and Eq. (<ref>) appear to be identical. However, this is not the case. The point of the definition of the natural orbitals is that they are expressed in a single-particle basis (only one coordinate 𝐫 as used more extensively in Apx. <ref>). Meanwhile, the eigenvectors used in Eq. (<ref>) have one coordinate 𝐫∈ℝ^3 for each electron in the system.  

We will note that certain renormalization schemes can generate a more efficient basis set than natural orbitals <cit.>. However, reducing the problem down to few enough orbitals that a polynomial time solver could be used would imply that the determination of that transformation is not discoverable in polynomial time since the general problem is known to be hard <cit.>. Thus, no universally efficient procedure should be expected.


For a given Hamiltonian operator H, the expectation value E (energy) is

    E=Tr(ρ H)≡∑_k⟨ k|ρ H|k⟩.

Replacing H by any other operator W gives the expectation value ⟨ W⟩. The index k is taken over any set of orbitals that is orthogonal and complete.


Throughout, we will only consider orthogonal basis states, which applies equally to k above, and the final result requires a local operator.

Typically what is done in practical computations is to take a truncated trace from the natural orbitals with the highest to the lowest weight. The result answer converges very quickly which can be seen from the using these functions in practice <cit.>.

For completeness, we define an excitation in the system using the above concepts.


Given a set of excitations spanning an interval of energy, the next excitation can be defined as follows. To find an excitation at eigenenergy E̅, take the set of all density matrices composed of eigenstates with E<E̅. Then find the states that are orthogonal to those states. The minimum energy will be the excitation up to degeneracy.


Two different excitations will have a different density matrices.



 §.§ Truncated density matrices


So far, we have discussed density matrices where the basis states used to describe the density matrix spans the entire Hilbert space. Let us now define a truncated density matrix where small occupations are set to zero.

Returning to Def. (<ref>), the definition can be modified to define a  truncated density matrix if the sum over k in Eq. (<ref>) is restricted to a value m less than the Hilbert space size, M.  We denote the truncated trace as

    Tr_m^(γ)(ρ^(α))=∑_i=1^m⟨Φ_i^(γ)|ρ^(α)|Φ_i^(γ)⟩

for orbitals from a set γ on the α excitation.


We use the term `truncated trace' because `partial trace' is usually associated with tracing over lattice sites and producing the partial density matrix.

There is an immediate consequence that the density matrix is now truncated, leading naturally to the definition of the truncation error.


The density matrix may be truncated to dimension m, known as the  bond dimension. The difference from the true value of the trace of ρ is known as the  truncation error, δ, which provides an estimate of the precision of the resulting expectation values from Def. <ref>. The full definition is then

    Tr(ρ^(α))-Tr_m^(α)(ρ^(α))=∑_i=1^Mρ_ii^(α)-∑_i=1^mρ_ii^(α)≡δ^(α)_m;γ

where α denotes a state that was used to construct ρ and γ is the basis over which the truncated trace was evaluated.


§ RELATIONSHIP BETWEEN TWO DENSITY MATRICES


The energy difference between excitations can be defined using the expectation values as the following.



For a given Hamiltonian H, the energy difference between two states α and β is

    Δ E_αβ = Tr((ρ^(α)-ρ^(β))H)

where Δ E_αβ=E_α-E_β. 


Two excitations α and β each have natural orbitals. These do not need to be the same set, nor orthogonal to each other. However, we restrict ourselves to orthogonal sets of natural orbitals so that they are related by a unitary transformation in every case.

The trace of any density matrix is invariant to the basis over which the trace is performed,

    Tr(ρ)=Tr(Uρ U^†)

this is often known as the cyclic property of the trace. However, a unitary cannot be assumed at the start, and we will see that the traces between the two density matrices must be satisfied to ensure this is true.


Natural orbitals of two different density matrices (assumed to be written in the same orbital basis with different occupational weights) are related by a unitary transformation if the trace is the same.

Two density matrices for two states α and β satisfy

    Tr(ρ^(α))-Tr(ρ^(β))=Δ N_e

and by the cyclic property of the trace from Eq. (<ref>), we can transform ρ^(α)→ Uρ^(α)U^† without changing the trace. If the number of particles between the two density matrices is the same, then this implies that the natural orbitals between the two states are related by a unitary matrix. For two states with different particle number, we would need to use an orthogonal transformation U^†→ O^-1 but the change in magnitude of the density matrix by the orthogonal transformation should be considered that we always consider natural orbitals with a normalized amplitude. Thus, one can enforce a normalization to again show that the natural orbitals are related by a unitary transformation.


















So long as the natural orbital sets for α and β are both complete in the space spanned by two density matrices (including unoccupied natural orbitals), then this proof holds. We restrict our consideration to density matrices with the same total trace for simplicity, but all results can be extended to the case where an orthogonal transformation is required.



 §.§ Local systems


The notion of locality can also apply to the eigensolutions and general operators <cit.>. In fact, is a central idea in quantum physics for all types of physically relevant ground-states.



  §.§.§ Local correlations

The area law is a statement of correlations for extremal eigenvalues of the full spectrum. The two conditions for correlation functions (operator O) depending on whether there is a gap in the eigenvalue spectrum (gapped) or not (gapless) <cit.>:

    ⟨ i|O|j⟩∼exp(-|i-j|/ξ)     [gapped]
    
    1/|i-j|^γ     [gapless]

for two arbitrary, real exponents ξ and γ. This condition holds up to a sign of the function. It is without loss of generality to more than two sites that this same definition still holds in the following.


For a proof of the correlation dichotomy, in some contexts known as Kohn's near-sightedness principle, we refer the interested reader to Refs. hastings2004locality,kohn1996density,prodan2005nearsightedness.


The natural orbitals are local, defined as a non-zero on a compact subset of 𝐫∈ℝ^3. 

Since they are derived from ⟨ĉ^†_iσĉ_jσ⟩ for each element of the density matrix. Upon diagonalization, a linear combination of these elements will be the result. Since a linear combination of local correlation functions are local themselves, the natural orbitals are also local.


  §.§.§ Local Hamiltonians


The specific property of the system that we want to study here is for local Hamiltonians ( i.e. those with finite extent). The basic assumption at a coarse level constrains the long-range behavior of the Hamiltonian is contained in the following definition.

A local Hamiltonian satisfies the following two properties in the thermodynamic limit.

    |i𝒜-jℬ|→∞lim⟨ i𝒜|H|jℬ⟩→0

where i and the set of coordinates 𝒜. 

If the interactions are local, then in one limit 

    |i𝒜-jℬ|→0lim⟨ i𝒜|H|jℬ⟩→ C_i𝒜

where C_i𝒜∈ℂ is finite and real. This definition holds whether the sites i𝒜 and jℬ represent single sites or clusters of sites, but it represents the ultra local limit where the Hamiltonian appears truly local.


 §.§ Relationship between natural orbital states
 



It is useful to explore the relationship between the natural orbitals of two different states and how the unitary that connects them can appear. There are two broad categories that the unitary can take and it is worth explicitly defining each.  After defining the two cases, we remark on some physical cases where these can be found.


A unitary matrix U relating the natural orbitals of two excitations is nearly the identity except for a sub-block over m states if the two states α and β both have small truncation error in the m most important orbitals to ρ^(α).

Denote the occupation values of the rotated density matrix U ρ^(β)U^† as (Tr^(γ)_m denotes a truncated trace in a basis γ for m of the most relevant orbitals)

    Tr^(α)_m(Uρ^(β) U^†)=∑_k=1^m_kk^(β)

and the unrotated value would be the same but without the bar applied. The bar denotes the truncated trace in the space of the most relevant orbitals for α. In general, this does not have to be the same evaluation as over the most important m orbitals for β.

The following statement is equivalent to having a small truncation error in the most relevant orbitals for a state α:

    Tr(ρ^(β))-Tr^(α)_m(ρ^(β))=δ^(β)_m;α!≅δ^(α)_m;α≪Tr(ρ^(α))

where it is clear that the truncation error is small. The trace is taken over the largest 1<m M orbitals relevant for the most important natural orbitals for α, {Φ^(α)}_m.  If we have δ^(β)_m;α≅δ^(α)_m;α,then the most important states for α give a small truncation error for β.

In order for the condition in Eq. (<ref>) to hold, the following must also be true

    ∑_k=1^mε_k^(α)≈∑_k=1^m_kk^(β)=∑_k,k'=1^m∑_ℓ=1^Mε_ℓ^(β)u_ℓ ku_ℓ k' ^*δ_kk'

following Lemma <ref> for an element of the unitary u_kℓ. We use the index k for the basis of natural orbitals for α. The indices ℓ and ℓ' for the natural orbitals for β.

There is a relationship between the unrotated coefficients and the β state and those for the truncation in the states relevant for α (indexed by k),

    ∑_k,k'=1^m∑_ℓ=1^Mε_ℓ^(β)u_ℓ ku_ℓ k' ^*δ_kk'≤∑_ℓ=1^mε_ℓ^(β)

where the sum ℓ over the diagonal elements of ρ^(β) is taken to be over the most relevant m natural orbitals for the β state. The equality is true when the truncation error is small for both α and β and only relevant orbitals for α are allowed, then those same orbitals must be relevant for β. Thus, the terms in the identity matrix appear strongly diagonal except for an m× m block for states of low truncation error.


States where the unitary is close to the identity in the irrelevant (M-m)×(M-m) block of U will be called  similar excitations.



Two excitations have a similar set of m natural orbitals if their density matrices have the relationship that a unitary U transforming one natural orbital set into another has approximately the decomposition U=W⊕ P for a unitary W of size m× m and an second unitary matrix P of size (M-m)×(M-m).


This definition is not meant to be exhaustive or tight for all possible similar states, but it is sufficient for the states here.

The opposite would be  dissimilar excitations.


Two excitations have a similar set of m natural orbitals but are not approximately of the form U=W⊕ P as defined from similar states.


This type of state is what would be encountered in the case of a symmetry protected state. Alternatively, two very widely separated centers of a potential v(x) in an eigenvalue problem would also be dissimilar.



 §.§ Truncation errors as a metric distance


The set of truncated density matrices should also be discussed in the context of a normed vector space, but the metric must be defined appropriately. In essence, we ask: how closely related are two truncated density matrices?


The relative truncation error between two states α and β will be defined as

    r^(γ)_m(α,β)=|δ^(α)_m;γ-δ^(β)_m;γ|

where it is very important to note that both truncation errors were evaluated in the same m-sized basis generated from excitation γ (either α or β or some other state).


In order for the truncated set of natural orbitals from one excitation to be mapped onto a vector space with respect to the natural orbitals of another excitation, there must be a definition of the metric distance. We can show that the relative truncation is a suitable metric.


The relative truncation error can be used as a practical metric distance between states in almost all cases of practical interest. 

The absolute value of the truncation error, Eq. (<ref>), satisfies all necessary properties of the metric distance in almost every case of relevance. The metric r^(α)_m for some truncation error δ^(α) satisfies

  * r^(γ)_m(α,α)=0
  * r^(γ)_m(α,β)≥0
  * r^(γ)_m(α,β)=r^(γ)_m(β,α)
  * r^(γ)_m(α,β)≤ r^(γ)_m(α,ζ)+r^(γ)_m(ζ,β) [Triangle ineq.]

which are the properties of a metric <cit.>. The triangle inequality is satisfied in the same way that the ℒ^1-norm is constructed normally.  Points 1 and 2 follow by definition. However, it is technically possible that r^(γ)_m(α,β) can be equal to zero as is especially evident when m=M as the truncation errors are both zero. Yet, it is generally the case that Point 2 holds for general states even when evaluating the truncation error out to numerical precision for 1≪ m M. This is why we call it a practical metric instead of simply a metric. The symmetry condition (3) is satisfied if the orbitals used to evaluate the truncation error is the same for both states. The triangle inequality (4) can be readily verified.


The main takeaway from this identification of a metric is that the amount to which the two states differ (with respect to a common set of states γ) from each other is communicated through the truncation error. When the description of a given state is accurate to δ^(α), then the states that will be most efficiently bundled are those also with low truncation error δ^(β).

The practical effect of this is that the truncation error is not only meaning the amount of information lost in a truncation to m orbitals, but the difference between truncation errors is also communicating the distance between the excitations, a highly remarkable feature especially in the context of entanglement renormalization which normally only assigns the truncation error as a means of uncertainty. The metric here implies it can also determine how well an excitation is described in a given basis.



 §.§ Energy differences as a metric distance


With regards to the identification of the metric for the truncation error, it is well-known that the truncation error implies an uncertainty in the energy, so it is often used as an error measure for entanglement renormalization methods <cit.>. Because of this relationship between the truncation error and uncertainty on expectation values, it is natural to ask if another quantity can also serve as a metric distance, a quantity that is more physical.

The most readily available quantity in many cases is the energy difference. We seek to establish if the differences in energies can be proven to be a metric as well.


Two density matrices (for states α and β) truncated to order m in some basis have a metric distance given by the absolute value of the energy difference, |Δ E^(m)_αβ;γ| with Δ E^(m)_αβ;γ=E^(m)_α;γ-E^(m)_β;γ, for local Hamiltonians expressed in local basis sets if not both of E_α and E_β are zero.

Consider the energy difference in a given basis written as

    Δ E_αβ;γ^(m)=.Tr_m^(γ)((ρ^(α)-ρ^(β))H)|_γ=α

where we select the γ=α basis here for clarity but can select any basis without loss of generality. An equivalent way to express this is to take the trace over all real-space positions but only expand the density matrices out to order m in the most relevant orbitals for α. The resulting energy difference in the most relevant natural orbitals for α is then

    Δ E_αβ;α^(m)=   ∑_i𝒜⟨ i𝒜|∑_k=1^m(ε^(α)_k|Φ^(α)_k⟩⟨Φ^(α)_k|.
       .-∑_ℓ=1^M∑_k'=1^mε^(β)_ℓ u_ℓ ku_ℓ k' ^*|Φ^(α)_k⟩⟨Φ^(α)_k'|)H|i𝒜⟩

where ρ^(β)=∑_ℓε_ℓ^(β)|Φ_ℓ^(β)⟩⟨Φ_ℓ^(β)| with ℓ indexing the complete space and with the unitary defined earlier in Eq. (<ref>).

 Let the definition

    _kk'^(β)=∑_ℓ=1^Mε^(β)_ℓ u_ℓ ku_ℓ k' ^*

hold in the following. A complete set of states can be inserted before H to give

    H|i𝒜⟩→∑_jℬ|jℬ⟩⟨ jℬ| H|i𝒜⟩

where the term ⟨ i𝒜|H|jℬ⟩ can be thought of as the Hamiltonian in real space. For local Hamiltonians, this term goes to a delta function for large differences between the sites i,j,𝒜,ℬ denoted as h_ij𝒜ℬ. 

In the local limit where (h_ij𝒜ℬ→ C_i𝒜δ_i𝒜,jℬ), Eq. (<ref>) becomes

    Δ E_αβ^(m)=   ∑_i𝒜⟨i𝒜|∑_k=1^m(ε^(α)_k|Φ^(α)_k⟩⟨Φ^(α)_k|.
       .-∑_k'=1^m^(β)_kk'|Φ^(α)_k⟩⟨Φ^(α)_k'|)|i𝒜⟩ C_i𝒜


At this point, the density matrix for β projected into the orbitals for α is not diagonal. If we impose that the orbital basis is local, as in Lemma <ref> for natural orbitals, then this implies that the as the difference between k and k' are large, then the overlap between functions must go to zero. Recall that we only consider orthogonal basis sets. This means that ∑_ℓ u_ℓ ku_ℓ k'^*=δ_kk' in the ultra-local limit. Consequently, the expression reduces to 

    Δ E_αβ^(m)=   ∑_i𝒜∑_k=1^mC_i𝒜(ε^(α)_k-^(β)_k)|⟨Φ^(α)_k|i𝒜⟩|^2.

where _k^(β) is the diagonal coefficients of ρ^(β) in the basis of the most relevant orbitals for α. 

With an absolute value sign, |Δ E^(m)_αβ;γ|, one can verify the axioms of a metric for 1≪ m M are practically satisfied as presented in Thm. <ref>. 

There can be an ambiguity in the definition of the metric if the states are zero but this is not typical of eigenvalue problems relevant here. So, we exclude this case.


Theorem <ref> made heavy use of the ultra-local limit where the Hamiltonian is effectively diagonal. This and several other features are discussed here as useful to derive the core statements, but understanding how things would change if the Hamiltonian is longer-range is worth performing. 

If the Hamiltonian contains longer range terms, then terms like ⟨ i𝒜|H|jℬ⟩ to some relative distance between i𝒜 and jℬ can be incorporated.


Many common basis functions decay exponentially. Notably, the Gaussian basis set decays exponentially away from the origin of the Gaussian. Thus, in the regime where matrix product states are best applied (local, gapped Hamitonians), the local approximation is a very good one for the excitations.

The ultra-local limit is not the most general form that the unitary can take here for arbitrary problems, but it is most applicable one for the tensor network case we find below where there is merely a permutation of the elements in the same basis. All of the excitations will be written in the same basis of entangled states and the density matrix eigenvalues of each is in an m× m basis ( i.e., effectively all of the occupational weights can be discovered in g same-size matrices which all are diagonal and in the same basis). If one prefers, then {Φ_k^(α)} and {Φ_k^(β)} belong to a common set γ and an m× m tensor with occupational weights can be identified for both states in that same basis.



 §.§ Large energy differences


We can now make a general statement about large energy differences.


Large energy differences imply larger differences between density matrices.

Note that

    |Δ E_αβ;γ^(m)|≤|Tr_m(|ρ^(α)-ρ^(β)|H)|≤Tr(|ρ^(α)-ρ^(β)|H)

where |ρ^(α)-ρ^(β)| is the element-wise absolute value implemented in the following way

    |ρ^(α)-ρ^(β)|=∑_k,ℓ=1^m|ρ̃_kℓ^(α)-ρ̃_kℓ^(β)||Φ^(γ)_k⟩⟨Φ^(γ)_ℓ|

where tildes denote that the density matrices are rotated into the γ basis. Effectively, one simply replaces round braces around the occupation values in Eq. (<ref>) by an absolute value. By comparison with Eq. (<ref>), then Eq. (<ref>) is satisfied. Further, recall that coefficients C_i𝒜 in Eq. (<ref>) set an effective energy scale. The maximum of which, C_max., can act as a normalization to give

    |Δ E_αβ;γ^(m)|/C_max.≤|Tr_m(|ρ^(α)-ρ^(β)|H)|/C_max.≤|ρ^(α)-ρ^(β)|_F

where|ρ^(α)-ρ^(β)|_F is the Froebenius norm and since coefficients of H can be either positive or negative. This is true by inspection of Eq. (<ref>) and |ρ^(α)-ρ^(β)|_F=∑_k,ℓ=1^m|ρ_kℓ^(α)-ρ_kℓ^(β)|. The sum over k and ℓ in the Froebenius norm can be truncated to m or kept in the full basis.

Summarily, when |Δ E_αβ;γ^(m)| is large, then the density matrices must be different as conveyed by the Froebenius norm. So, a large energy difference implies a large difference between the states α and β.


The result is valid in any basis if the operator is local. The use of a finite number of states m for the orbital basis makes this result useful in a truncated space.

So, density matrices with large energy differences necessarily have large differences. The opposite is not so clearly defined. If the energies are low, they can either be a similar state, in which the most relevant basis for one excitation is relevant for another state. Alternatively, if the state is dissimilar, then the basis relevant for both excitations is very different, as discussed earlier.

There are several useful consequences of the previous derivation of the main result. We will discuss them here before moving towards the solution with partial density matrices.


The results of Thm. <ref> generalize to all nearby (Δ E_αβ≈0) states with low truncation error.

Reconsider the form of Eq. (<ref>) but now for a chain of energy differences between several states of a single symmetry sector. The first energy difference satisfies Eq. (<ref>) as does the next set of two excitations. Therefore, as a transitive property, all wavefunctions with small energy differences share a high degree of overlap in their truncated density matrices.


Taken together, the results of this section show that partial density matrices with low energy can be bundled efficiently ( i.e., a basis m can give a small truncation error for both states) together. Each excitation added to the bundle comes at either no additional cost in terms of the number of natural orbitals required for a small truncation error if the states are similar. For small energy differences, there is a possibility that the small energy states cost no more than making the lowest energy state. When bundling two states in the bulk of the eigenvalue spectrum together, one must only pay the cost of one of the states, the second comes at a small cost if similar.

When adding a dissimilar state to the bundle, a sufficient number of natural orbitals must be added to the orbital basis to give a low truncation error for that dissimilar state. Once the new dissimilar state is added to the bundle, a new set of nearby similar states can be bundled at small cost.



 §.§ Application to quantum chemistry


The primary focus of this paper is on locally entangled systems. This is not the general case in quantum chemistry where models have longer range interaction. In this case, the entanglement of the states is known to be much larger than for the local models as can be seen from direct computation. In particular, we note that the extension of models to two dimensions can cause an exponential increase in the entanglement with the width of the system <cit.>.

This would cast doubt on whether the ultra-local limit will apply in a quantum chemistry system. It turns out that by use of the singular value decomposition (SVD), the basis states of the model can be written such that the basis functions between two states are identical for two different systems. For example, writing two excitations and then cutting the system with the SVD at the same bond in each system will give a form like Ψ=UDV^† which has grouped basis functions to the left of the cut (contained in U) and the other basis functions for the right of the cut (contained in V^†). Thus, the same basis functions can be assigned for both U and V^† for both states. Between the two states, the basis functions may be different in practice, but we can always find a unitary transformation that makes the states match between the U (V^†) for one state and the other.

However, the two D matrices are still diagonal. This is exactly the ultra-local limit as we applied it to the density matrix (but not the Hamiltonian from Def. <ref>) because the only difference is that the two states have different occupation numbers inside of D. Thus, only a permutation is allowed from the basis functions to rearrange the occupation numbers of the first state to the second. Thus, the ultra-local limit, in the construction of the SVD, is valid even between two different states. We merely motivated the limit by applying an understanding of the real-space orbitals previously.

Thus, quantum chemistry can fit into this hierarchy established here. The computational bottleneck appears in retaining enough states in the SVD in the quantum chemistry system. Because entanglement is larger, there must be more states retained, thus the method is less efficient.



§ DENSITY MATRIX RENORMALIZATION

Up to now, all considerations have been for the full wavefunction and full density matrix and n-body reduced density matrices. The statement for the full wavefunction is not completely useful for the solution in a tensor network decomposition where partial density matrices on either the left or right partition of a given system (partitioned in the sense of partitioning a graph) is the relevant quantity for the eventual solution. 


A system partitioned into two sets can identify basis functions for each set. The basis functions may be used to project the full density matrix into a reduced site representation. The density matrix on one half of the system is equivalent to the full density matrix but traced out on the other half. No matter how the system is partitioned, the entanglement expressed by one partial density matrix must be equal to the complementary density matrix's entanglement.


We do not rule out that a null set can be used for one of the two partitions, but this would imply that the partial density matrix is equivalent to the full density matrix. The sites do not need to be contiguous, but they are chosen to be such here.



 §.§ Matrix product states


An entanglement renormalization algorithm explicitly calculates the components of a density matrix partitioned between two parts of a system with the use of a singular value decomposition (see Ref. bakerCJP21,*baker2019m,dmrjulia1 for an explicit derivation). The following definitions cover the matrix product state and make use of the above theorems.


A wavefunction with degrees of freedom σ_i on each site i is written as

    |ψ⟩=∑_{σ_i}c_σ_1σ_2σ_3σ_4…|σ_1σ_2σ_3σ_4…⟩

for some probability amplitudes c_σ_1σ_2σ_3σ_4…∈ℂ. By performing a series of reshapes and SVDs <cit.>, Eq. (<ref>) can be decomposed into a series of tensors <cit.>
    |ψ⟩=∑_{σ_i}, {a_i}
    A^σ_1_a_1A^σ_2_a_1a_2… D^σ_i_a_i-1a_i… 
    B^σ_N_a_N-1|σ_1…σ_N⟩

for a number of sites, N. Raised and lowered indices are for ease of viewing as there is no notion of covariant or contravariant indices <cit.>. Raised and lowered indices are used here to signify the partitioning of the lattice via a reshaping operation <cit.>. The index introduced by the SVD (a_i) is known as a link index and its dimension is called a  bond dimension. The orthogonality center D contains the weight of the basis sets and can be gauged to any site or bond <cit.>. All tensors to the left of the center of orthogonality (A) are  left-normalized. Similarly, all tensors to the right of the orthogonality center (B) are said to be  right-normalized. Contraction of any left- or right-normalized matrix with itself leaving only two indices corresponding to a_i uncontracted yields an identity matrix.


Eq. (<ref>) can be represented with Penrose's graphical notation as in Fig. <ref><cit.>. Vertical lines on the MPS correspond to the σ indices and a are the horizontal indices.

Because of how the density matrices are truncated in the MPS through the D matrix retaining the largest values, we can remark that the basis states that are kept represent the most entangled basis functions between the two partitions. Thus, the use of natural orbitals in the previous derivations and theorems now becomes a set of the most entangled basis functions between the two partitions.



 §.§ Matrix product states represent the ground-state faithfully


In a seminal work in Ref. verstraete2006matrix, it was demonstrated that the MPS can represent the true ground state with a low error. The rigorous result is established with a bound on the Renyi entropy in relation to the truncation error of the MPS. The argument demonstrates that locally entangled models are efficiently described by the MPS ansatz and there is a generalization to the MPS in higher dimensions and connections to the area law of entanglement.

The techniques used in this paper veer closer to those used in quantum chemistry, making use of the MPS's relationship to the full density matrix. We recast the results of Ref. verstraete2006matrix into the tools used here in case it is useful.

Consider a wavefunction of N sites. Reshaping the N degrees of freedom into two groups, a left group and a right group allows us to take a singular value decomposition of the form (we drop the dagger from V when writing it in terms of tensor components for clarity)

    ψ_σ_1…σ_N=U^σ_1…σ_j_a_j-1D_a_j-1a_jV_a_j^σ_j+1…σ_N

where σ_i represents the degrees of freedom locally on each site and the wavefunction was decomposed on the jth bond. Raised and lowered indices do not mean anything and are only used for clarity.

The decomposition according to the SVD gives the elements necessary to construct the partial density matrix for the left

    ρ_L=UD^2U^†

and right

    ρ_R=VD^2V^†

of the system <cit.>. The occupations of the natural orbitals are contains in the D^2 matrices and are known to decay rapidly. Note that if we had the expectation value of the Hamiltonian, H,

    E=Tr(ρ H)=∑_iλ_ih_ii

where the basis of the natural orbitals was used to make the density matrix diagonal. If the orbital occupations of the density matrix are ρ_ii=λ_i and ordered as 1≥λ_1≥λ_2≥λ_3≥λ_4…≥λ_N≥0, and decay rapidly, which is an assumption we use throughout.

Because the occupation values of the density matrix decay rapidly, the entanglement, S=-Tr(ρlnρ) is well approximated. Thus, if the entanglement is low, then the truncation error is also low <cit.>.

This argument we use here leaves out the extensions to the area law of entanglement, but we find these methods useful for the extension to the bundled MPS case.



 §.§ Bundled matrix product states


The bundled MPS represents several density matrices considered together but independently. We note a key distinction that when writing a wavefunction for two excitations Ψ=(ψ_1+ψ_2)/√(2) that the states individually would form a pure density matrix but that the combination is mixed. The density matrices for each state in the bundle would still be pure. Since the wavefunctions are considered independent from each other is the key difference between the bundled and ensemble density matrices.  The most common form will be a set of mixed density matrices (with or without normalization of the trace to a value of 1).



A ensemble of excitations can be added to the MPS with the additional of another index ξ onto the orthogonality center.

    |Ψ⟩=∑_{σ_i},
    {a_i}, ξ
    A^σ_1_a_1A^σ_2_a_1a_2… D^σ_iξ_a_i-1a_i…
    B^σ_N_a_N-1|σ_1…σ_N⟩|ξ⟩

where ψ_ξ represents g∈ℕ excitations indexed by ξ. This is also referred to as a bundled MPS.


 §.§ Construction of the bundled matrix product state


An explicit construction of the exact bundled MPS was demonstrated in Ref. baker2024direct from a group of MPSs. Note that that construction in Ref. baker2024direct over-determines the form of the left- and right-normalized matrices. The ability to compress the bundled MPS is the main point illuminated by the results in this paper ( i.e., many states are similar enough to reduce the size of the left- and right-normalized tensors).

An alternative derivation to Ref. baker2024direct can be found from the traditional form of a set of excitations as an d^N× g matrix for N sites, size of the physical index d, and number of excitations g, giving a complex coefficient, c_σ_1σ_2…σ_N^ξ where σ_i is defined as in the main text to index the physical degrees of freedom on a site i and ξ indexes the excitations.

The task to separate the indices follows the same procedure as the derivation of the MPS from a single wavefunction <cit.>. The first step is to reshape the matrix to isolate σ_1 from the rest of the indices and perform a singular value decomposition.

    |Ψ⟩=∑_σ_1…σ_N
     a_1ξA^σ_1_a_1c^ξσ_2…σ_N_a_1|σ_1…σ_N⟩|ξ⟩

Continuing with more SVDs and reshapes, we discover more terms in the bundled MPS until recovering Eq. (<ref>). 

The importance of this exercise is to discover the size of the bond dimension at the ends of the bundled MPS.  Recall that for the single MPS that the maximum size of the bond dimension in the exact case (without truncation) scales as

    m_max=min(∏_x=1^id_x,∏_x=i+1^N_sd_x)

for bond i and physical index of size d_x on bond x. This is because the SVD is chosen with a convention that the size of the D matrix is the minimum of the two dimensions of the input matrix <cit.>. 

A similar expression can be derived for the bundled MPS. The only difference is that when gauged to a site j, the excitation index can be thought of as multiplying the bond dimension for the purposes of this size counting argument. Thus, redefining Eq. (<ref>) with d_j→ g· d_j will be sufficient.

The main point of this exercise is to show that when the orthogonality center is gauged to the first or last site in the system, one half of the system has the same size bond dimensions as the single MPS. The other half of the system has a larger bond dimension dependent on g.

This procedure is unfeasible beyond small systems since the exact wavefunctions can not be determined efficiently by exact diagonalization from the exponential size of the Hamiltonian operator with increasing number of sites <cit.>.



 §.§ Overlap matrices of system partitions


The overlap matrices representing the left and right partitions of the system can be constructed from the left- and right-normalized tensors of the MPS.


The wavefunction ψ may be decomposed with an SVD as (with the matrix D given on a bond, not a site)

    |ψ⟩=U^σ_1σ_2…σ_i_a_i-1D_a_i-1a_iV^σ_i+1…σ_N_a_i|σ_1…σ_N⟩

The contraction of the wavefunction onto another state ϕ is given by

    ⟨ϕ|ψ⟩=∑_a_i-1,a_i,a'_i-1,a'_iρ^(L)_a_i-1a_i-1'D_a_i-1a_iD̃_a_i-1'a_i'ρ^(R)_a_ia_i'

and with

    ρ^(L)_a_i-1a_i-1'=∑_{σ_i}Ũ^σ_1σ_2…σ_i_a_i-1'U^σ_1σ_2…σ_i_a_i-1
    ρ^(R)_a_ia_i'=∑_{σ_i}V^σ_i+1…σ_N_a_iṼ^σ_i+1…σ_N_a_i'

with both of those terms represented in Fig. <ref>. The matrices D can be truncated just as the density matrix was to m≤ M just as in a truncated SVD. This will define a  truncated partial density matrix.


The following minor statement is presented as a lemma because it directly follows from the additional definition of the partial density matrix and Thm. <ref>.


A set of relevant singular vectors for a partial density matrix of an excitation can differ at most from the most relevant singular vectors of another excitation related to their energy difference similar to Thm. <ref>'s statement for natural orbitals.

The partial density matrices are presented graphically in Fig. <ref> between the excitations.  Note that immediately, one can apply the results from Thm. <ref> onto the partial density matrices.  For those states on the left of the orthogonality center, the density matrix is written in the states for the left and right basis functions

    ρ=
     U  D^2 U^†   (left)
    
     V  D^2 V^†   (right)

and the indices on each tensor are skipped here for the sake of brevity.

In order to use Eq. (<ref>) with the density matrix, the lattice sites on the whole system must be partitioned into a left (L) and right (R) group ( i.e., either (i𝒜)_L or (i𝒜)_R). The sum can then be separated into each group, and the density matrix must be represented in a basis natural to each group ( i.e. with the basis functions for the left or right as in Def. <ref>.

    Δ E^(m)_αβ;γ   =∑_j∈{(i𝒜)_L}C_j⟨ j|U_γ(D_α;γ^2- D_β;γ^2 ) U_γ^†|j⟩
       =∑_j∈{(i𝒜)_R}C_j⟨ j|V_γ(D_α;γ^2 - D_β;γ^2)V_γ^†|j⟩

where the subscript denotes which excitation the component of the SVD belongs to for some chosen basis γ of which m states were kept. Just as in Eq. (<ref>), the Hamiltonian must be local for this expression to be valid. In the tensor network formalism, this corresponds to the MPO containing only local terms. This is because the states k from Eq. (<ref>) can also index the basis vectors of the SVD used to compose U, D, and V. This proof applies to any bond in the MPS, so there is no issue of how the MPS is gauged.


Note that the D matrices for either the α or β excitation are guaranteed to be diagonal by construction and are written in the same basis of entangled states. This is exactly the ultra-local limit from Def. <ref> that was used in Thm. <ref>. In this case, the orbitals {Φ^(α)_k} and {Φ_k^(β)} are the same basis. The degree of freedom that allows one to write a higher excitation is the changing of the magnitude of the occupancies ε_k for either state. Thus, the general prescription defined for the general case of natural orbitals in the quantum chemistry context is precisely the relevant case for the tensor network.

For smaller energy differences, the most relevant singular vectors will be more common when the states are similar. This means that the primary results from Thm. <ref> can be equally applied to the singular vectors here implying a maximum bound on the number of differences in D. The singular vectors are also orthogonal in this case, so this is different from Thm. <ref> in that the overlap of the singular vectors is not needed.

The reduction here applies equally well for pure or mixed states, although in a given computation the states computed are a mixed state representation of a given eigenvalue. Combined with the previous results, the main statement here applies to all local Hamiltonians.

We are then ready to state a central theorem to the paper.

A bundle of states in a tensor network containing an energy interval [E_α,E_β] (Δ E_αβ≈0) and with δ^(α)≪1 and δ^(β)≪1 will not require a larger bond dimension m to add a similar state δ^(γ)≪1 with Δ E_αγ≈0. It will cost at most m' more states (the number required to describe γ with δ^(γ) of some magnitude) to add a dissimilar state of any energy difference.

This immediately follows from the results in Sec. <ref> since the energy of the new state on the interval will have a small energy difference with the other states on the compact interval. The summary statement is that the bundled MPS is not necessarily more expensive to solve for if a suitable initial state is solved and then low energy excitations are found.



 §.§ Discussion of states with different symmetries



How many excitations are similar is highly model dependent. While it is true that the next excitation does not need to be similar, it is true that the set of all states on a compact interval contains similar states in most models. So, there should be an expected reduction in the size of the bond dimension.
 
 Consider two classical MPSs (m=1), one with a state with 1 fermion and another with 10 fermions. The overlap between these two states will definitely zero because of the quantum number symmetries. As the bond dimension of the MPS is grown (as in a DMRG computation) and each tensor contains more blocks each with different symmetries, it is not guaranteed that the symmetry sectors will have much overlap upon contraction of the full network. From this argument, it can be reasoned that based on symmetries alone that the singular values will not have much overlap between the two excitations. This is even true for a state with 1 fermion and another with 2 fermions with all fermions on the first site and if the net fluxes of the quantum number symmetries are assigned on the last site.

Conservatively, we can state that the addition of a new state into the bundled MPS that does not have a symmetry represented will require another m states to be added onto the link index. So, the scaling appears as O(mN_sym.); however, this a generous upper bound because of similarities between some symmetry sectors.

Take for example the ensemble containing 1 spin-up fermion and another excitation for 1 fermion with a spin-down fermion for a spin-symmetric Hamiltonian ( i.e., no magnetic field or other effect applied). These two symmetry sectors in the Hamiltonian are identical blocks by definition. The wavefunctions, however, while not having the same symmetry have many elements in the Hamiltonian in common and have a conjugate set of singular vectors in the other.

Another example would be for a chain of entirely up spins and another state of entirely down spins. The quantum number symmetries on each link index show that the spaces in which the tensors exist is completely orthogonal to each other in general.



§ EXAMPLES


To illustrate some of the concepts presented in the previous sections, we focus on three models and determine some relevant quantities. All models are solved with the DMRjulia library.



 §.§ Transverse field Ising model


We consider the case of the transverse field Ising model, defined as

    H=∑_iσ_i^z·σ_i+1^z+h_xσ^x_i

where

    σ^x=([ 0 1; 1 0 ])  and  σ^z=([  1  0;  0 -1 ])

with subscripts indicating the matrix belongs to a Pauli string of Kronecker products O_i=I⊗ I⊗…⊗ O⊗…⊗ I with O in the ith position. The model exhibits a phase transitions when h_x=1, which creates a gapless spectrum <cit.>. We note that we focus on the properties of the bundled MPS in this paper and not the algorithm used to solve it, which in the most general case should be considered to be undecideable–in the sense of the halting problem–to find the spectral gap of the Hamiltonian <cit.>, meaning that finding the right method to solve the bundled MPS is an open challenge for arbitrary systems (including three-dimensional models).

One way to visualize the relevant states in the system for different energies is to form the overlap matrix, Γ from Fig. <ref> and multiply the D matrices contracted onto the link index, giving

    Γ_ij = Γ_a_i'a_i=∑_a_i-1,a_i,a'_i-1,a'_iρ^(L)_a_i-1a_i-1'D_a_i-1a_iD̃_a_i-1'a_i'

This ensures that no only will we be able to see the overlap between two of the basis states used in two different bundles, but we will also see how relevant each are for the overall answer. All plots are shown with the logarithm of the absolute value of each element of the overlap matrix Γ. Matrices are truncated so that no singular values of zero are represented. 

Brighter colours indicate that the matrix element is important and should not be truncated. Darker colours indicate that the matrix elements is small. Truncating an entire row or column is possible in the normal evaluation of the reduced bond dimension of a system. Thus, we look in the figures for a row or column where (nearly) all of the elements are a darker colour. This means that the row or column can be removed without losing precision in the bundled MPS.

The results in Fig. <ref> tell the same story as the theorems presented previously. The larger the energy difference between states, the more states in this measure have a higher weight (lighter color). 

The general trend is that the larger the energy difference, the fewer rows and columns that can be truncated.

Results for the critical model in Fig. <ref> are largely the same but the bond dimension is larger because of the gapless eigenspectrum of the model.



 §.§ Heisenberg model


We can perform the same analysis with an XXZ model (Δ=1) of the form

    H=∑_i𝐒_i·𝐒_i+1

where 𝐒=⟨ S^x,S^y,S^z⟩=1/2⟨σ^x,σ^y,σ^z⟩ and the commutator [σ^x,σ^z]=-2iσ^y<cit.>.



Figure <ref> shows the weighted overlap matrix. This time, groups of 10 eigenstates are shown. The Heisenberg model has a well-known symmetry structure that creates a blockier appearance for the graph. However, the main message of the analysis remains the same: larger energy differences have more high-weight states. Bundling together energy states that are close together have many states in common and thus the bond dimension can be truncated.



§ CONCLUSION


We introduced the bundled density matrix, a set of density matrices that are independent but written in a common basis. We showed that the difference in the truncation errors constitute a practical metric to determine a notion of distance between the density matrices. It was shown that for local system, the energy difference in that basis was also a practical metric to describe a notion of distance between density matrices in the bundle.

The larger the energy difference between the density matrices, the larger the difference in the density matrices are. One might expect that as one bundles excitations that are further into the bulk that the volume law entanglement would dominate and therefore drive the bond dimension higher. What these results suggest, in effect, is that this is only guaranteed for large energy differences. For small energy differences, the density matrices can either be similar, where the m most relevant states constitute a good basis for low-energy excitations. Or, the states can be dissimilar, such as in the case of a different symmetry sector, where the truncation error is large.

We extended these ideas to the bundled matrix product state where the results demonstrate that the bond dimension is lower for similar states. Bundled matrix product states therefore may support low energy excitations in the same symmetry sector without a costly increase in the bond dimension. To add a dissimilar state, one must pay the cost of the matrix product state representation of that state. Large energy differences in bundled matrix product states have effectively no degrees of freedom overlapping and are not very different from two separate MPSs.

Definitions for similar states and the exact applicability of the local limit used in this paper can certainly be expanded and modified, and we encourage the community to work with the concept of the bundled density matrix in other contexts including quantum information.



§ ACKNOWLEDGEMENTS


This research was undertaken, in part, thanks to funding from the Canada Research Chairs Program.

The Chair position in the area of Quantum Computing for Modelling of Molecules and Materials is hosted by the Departments of Physics & Astronomy and of Chemistry at the University of Victoria.

N.S. acknowledges the NSERC CREATE in Quantum Computing Program (Grant Number 543245).

This work is supported by a start-up grant from the Faculty of Science at the University of Victoria. 

This work has been supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC) under grants RGPIN-2023-05510 and DGECR-2023-00026.



§ TWO EQUIVALENT FORMS TO DETERMINE THE DENSITY MATRIX


The definition of the one-body reduced density matrix is 

    ρ(𝐱,𝐱')=   ∫…∫ψ^*(𝐱,𝐫_2,…,𝐫_N_e)
         ×ψ(𝐱',𝐫_2,…,𝐫_N_e)d𝐫_2… d𝐫_N_e

which works well for applications in quantum chemistry. In the graphical notation used in the main section of this paper, the vertical physical indices would have to be kept track of as separate but all other physical indices would be contracted. 

The case of the tensor network has an immediately large computational cost with using this form. The tensor network form would require that we isolate physical indices and perform a contraction that winds up being exponential for large |i-j|. Neither of these is desirable for the tensor network because it was specifically formulated not incur the exponential cost of the full quantum problem.

Another definition of the density matrix with elements

    ρ_ij=⟨ψ| c^†_iσ c_jσ|ψ⟩

for fermions <cit.> is often cited. This form is equivalent to the above form as we now demonstrate.

Consider the definition of the wavefunction as a superposition in some position. We denote the local degrees of freedom on each cite i by σ_i. The general decomposition of the wavefunction becomes

    |ψ⟩=∑_σ_1…σ_Nw_σ_1…σ_N|σ_1…σ_N⟩

for N sites. The application of the operator c_j gives

    c_j|ψ⟩=∑_σ_1…σ_Nw_σ_1…σ_N|σ_1…σ̅_j…σ_N⟩

where σ̅_j is the modified value of σ_j upon evaluation of the operator. In the case of a spin-half system, a value of ↑ would go to ↓.

A similar form can be derived for c^†_i on the dual vectors

    ⟨ψ| c^†_i=∑_σ_1…σ_N⟨σ_1…σ̅_i…σ_N|w_σ_1…σ_N^*

Neither this sum nor the dual expression in Eq. (<ref>) have the same number of terms as in Eq. (<ref>) because the annihilation operator creates some zero terms from the original expression. The final expectation value becomes

    ⟨ψ| c^†_i c_j|ψ⟩=∑_σ_iσ_j⟨σ̅_iσ_j|W_σ_iσ_j|σ_iσ̅_j⟩.

where

    W_σ_iσ_j=∑_{σ_1…σ_N}\{σ_iσ_j}w^*_σ_1…σ_Nw_σ_1…σ_N

where the summation is over all variables but not σ_i or σ_j.
Returning now to Eq. (<ref>), it can be seen that the replacement of the 𝐱 and 𝐱' terms would be represented with the ket form as

    ρ(𝐱,𝐱')⇒   ρ(σ_𝐱,σ_𝐱')discretize=ρ(σ_i,σ_j)
    
    =   ∑_{σ_1…σ_N}\{σ_iσ_j}⟨σ_1…σ_N|W_σ_iσ_j|σ_1…σ_N⟩

and the primed index on i or j is a reinterpretation of the continuous, real-space position variable 𝐱 into the lattice variable for some site on the discrete lattice space. Eq. (<ref>) then reduces to Eq. (<ref>). Thus, the two forms of the density matrix are the same.