Efficient and Flexible Approach to Simulate Low-Dimensional Quantum Lattice Models with Large Local Hilbert Spaces

Quantum lattice models with large local Hilbert spaces emerge across various fields in quantum many-body physics. Problems such as the interplay between fermions and phonons, the BCS-BEC crossover of interacting bosons, or decoherence in quantum simulators have been extensively studied both theoretically and experimentally. In recent years, tensor network methods have become one of the most successful tools to treat lattice systems numerically. Nevertheless, systems with large local Hilbert spaces remain challenging. Here, we introduce a mapping that allows to construct artificial $U(1)$ symmetries for any type of lattice model. Exploiting the generated symmetries, numerical expenses that are related to the local degrees of freedom decrease significantly. This allows for an efficient treatment of systems with large local dimensions. Further exploring this mapping, we reveal an intimate connection between the Schmidt values of the corresponding matrix-product-state representation and the single-site reduced density matrix. Our findings motivate an intuitive physical picture of the truncations occurring in typical algorithms and we give bounds on the numerical complexity in comparison to standard methods that do not exploit such artificial symmetries. We demonstrate this new mapping, provide an implementation recipe for an existing code, and perform example calculations for the Holstein model at half filling. We studied systems with a very large number of lattice sites up to $L=501$ while accounting for $N_{\rm ph}=63$ phonons per site with high precision in the CDW phase.


Introduction
Large local Hilbert spaces appear in various kinds of problems in quantum many-body physics. Prominent examples arise in the field of ultra-cold quantum gases. Systems such as interacting bosons in a one-dimensional lattice [1,2] or trapped ion quantum simulators [3][4][5] have been studied extensively, fertilizing a rapid theoretical and experimental progress. Another typical problem featuring large local Hilbert spaces is the interplay between lattice fermions and phonons. For instance, the formation and stability of (Bi-)Polarons is a central problem and considerable effort has been taken for its investigation [6][7][8][9][10][11][12][13][14][15]. A broad class of different methods such as quantum Monte Carlo [16,17], density-functional theory [18], density-matrix embedding theory [19], or dynamical mean-field theory [20][21][22] has been explored to study its various aspects. Evidently, the task to numerically describe such low-dimensional, strongly correlated quantum systems has been subject to a vast development. In particular, the capabilities of tensor-network methods have improved a lot in the past two decades. Here, matrix-product states (MPSs) have become the fundament for flexible, numerically unbiased and in principle exact methods allowing for the study of not only ground-state properties but also of out-of-equilibrium dynamics of quantum many-body systems [23][24][25][26][27][28][29][30][31][32][33].
In (time-dependent) density-matrix renormalization group (DMRG) methods, the computationally limiting factor is the bond dimension of the tensors when performing tensor contractions [26-28, 30, 31, 33]. For instance, using MPS, one is mostly concerned with matrix-matrix contractions, which scale with the third power of the dimensions of the involved matrices. However, these operations can be rendered cheaper if the system under consideration conserves global symmetries. Being able to exploit (non-)abelean symmetries is an important feature of tensor networks in general [34][35][36][37], as a large bond dimension is related to the amount of entanglement and decay of correlation functions [29,38,39]. Aiming to describe strongly correlated systems, large bond dimensions can be required and thereby exploiting as many symmetries of the system as possible is highly desired.
Another important contribution to the numerical expenses of MPS algorithms is the dimension of the local Hilbert spaces H j . For instance, when considering systems with large spin or bosonic degrees of freedom, a local dimension dim H j ≡ d j ∼ O(10) − O(100) can yield drastic restrictions on the maximum possible bond dimensions as typical contractions usually scale with d 2 j or even d 3 j . In order to overcome such restrictions, approaches such as the pseudo site (PS) and the local-basis optimization (LBO) method [6,[40][41][42] were developed. These methods have proven to be successful tools for treating fermion-phonon couplings in the Holstein model, even out of equilibrium and at finite temperature [6,10,13,14].
In this paper, we introduce an alternative approach to simulate systems with large local Hilbert spaces efficiently and in a flexible framework. In order to treat these kinds of systems efficiently with MPS, we exploit the fact that global U (1) symmetries reduce effective local block dimensions drastically [34][35][36]. The starting point of our method is a thermofield doubling of the many-body Hilbert space, which is an established procedure in finite-temperature DMRG [43][44][45]. Then, introducing a new representation for operators in a particular subspace of the doubled Hilbert space allows us to show that global operators breaking U (1) symmetries can be identified with projected purified operators that conserve the corresponding symmetries 1 . Thereby, challenging general lattice systems breaking global U (1) symmetries with d j ∼ O(10) − O(100) can always be mapped into numerically more feasible systems. Importantly, this mapping requires only minor changes in existing codes and is completely general.
The paper is organized as follows: At first, we present the relevant aspects of our approach in Sec. 2 in a non-technical fashion and provide an implementation recipe that captures the changes in actual codes. In Sec. 3, we introduce the projected purification in great detail and show how to construct corresponding operators. In Secs. 4 and 5, we present the representation in terms of MPS. This includes the observation in Sec. 5.1 that connects the singular values of auxiliary bond indices with the diagonal elements of the single-site reduced density matrix, yielding an intuitive physical criterion to decide when the projected purification provides numerical benefits. We also discuss the numerical complexity of this mapping in actual calculations in Sec. 5.2. We close our discussion in Sec. 6 with an exemplary application of our mapping to physical Hamiltonians, i.e., the Holstein model and a Hubbard model with pair creation and annihilation and present some numerical results for the former.

General Concept and Implementation Recipe
The general idea of our mapping is to exploit global U (1)-symmetries, where the system under consideration does not conserve them in the first place. In the tensor-network framework, states can be constructed so that they transform under a global symmetry, i.e., they are eigenstates of the corresponding symmetry generator. Let us consider a system with a global particle number operatorN . Eigenstates |N ofN are labeled by their irreducible representations N and (ignoring degeneracies) any state can be decomposed in terms of these eigenstates doubling doubling projection projection Figure 1: Mapping of local operatorsb j andb † j acting on H into projected purified local operatorsb P ;jβ † B;j andb † P ;jβB;j acting on P. This transformation is the central, necessary modification for existing codes in order to use our method. Now, we can perform a doubling of the original Hilbert space and construct states of the form where we introduced labels P, B to distinguish the different Hilbert spaces. We restrict the allowed coefficients N such that each state |N can be mapped uniquely to a state |N P ⊗ |N 0 − N B with a properly chosen N 0 : The transformed wavefunctions are eigenstates of the new, global symmetryN P +N B with eigenvalue N 0 and can therefore be represented efficiently by symmetric MPS. The subspace P spanned by all states |N P ⊗ |N 0 − N B has the same dimension as the original Hilbert space so that no additional complexity is generated with this new representation. Importantly, the same considerations can be applied to the local degrees of freedom, constituting the many-body Hilbert space. Guided by this idea we will show in the following sections that there is a simple prescription to transform operators so that they are acting in P only. Using balancing operatorsβ ( †) B;j (which are introduced in Eqs. (17) and (18)), global operatorsÔ that break the global U (1) symmetry generated byN can be mapped into operators conserving the global U (1) symmetry generated byN P +N B . This is achieved by replacing ladder operatorsb The detailed mapping, containing also the intermediate step of doubling the Hilbert space, is shown in Fig. 1.
Recapitulating this short description of the general ideas of our mapping it should be noted that the states mapped to P are pure states in P but describe mixed states with respect to the orthogonal decomposition of H in terms of the eigenstates ofN . This is in close reminiscence to the purification procedure [43][44][45] that is commonly used to represent mixed states with respect to H. However, there is also an important difference: Restricting the allowed states by a projection into the subspace P of the doubled Hilbert space, the complexity of the state's representation is conserved, i.e., our mapping does not add additional degrees of freedom to the problem under consideration.
Next, we provide a short recipe, for how to implement the previously described projected purified DMRG (ppDMRG) for ground-state searches and time-evolution methods, including prerequirements. Note that this recipe is particularly short, because the necessary changes are small.
Prerequirements In order to incorporate ppDMRG into an existing framework, it is necessary that the framework can handle Hamiltonians with more than nearest-neighbor interactions.
Necessary changes The existing set of local operators needs to be extended with balancing operators that act on the bath sites, as introduced in Eqs. (17) and (18). In particular, for every species of local creation and annihilation operators corresponding balancing operators are needed when changing a global U (1) quantum number. Those operators shall only have zero and one as elements and always commute with every other operator. Additionally, for each species of creation-and annihilation operatorsb ( †) j , a parity-operatorPb j = e iπb † jbj might be useful. A scenario in which the action ofPb j is necessary is discussed in Sec. 6.
Usage Following these changes, all existing tools can be used as usual, but with a doubled system size where physical and bath sites alternate, which is a common technique in finite-temperature DMRG. Hence, local observables are now evaluated via two neighboring operators. However, care must be taken that the MPS represents states in P, i.e., the L local gauge constraints defined in Eq. (16) have to be fulfilled. Fortunately, since projected purified operators manifestly act on P only, it suffices to ensure that the initial state of any algorithm is in P. For instance, using the previous conventions, an initial state for a ground-state search is given by the product state Clearly, for typical ground-state calculations this state is a bad initial guess. However, it can be used as a starting point to create more suitable initial guess states by applying sequences of projected purified operators.

General Models and Bath Sites
We consider a lattice system of L ∈ N degrees of freedom, each of which being described within a Hilbert space H σ of local dimension σ ∈ N spanning the system's overall tensor-product Hilbert space H = H ⊗L σ . A state |ψ ∈ H can be expressed in terms of all local degrees of freedom |σ 1 · · · σ L ∈ H: with, in general, complex coefficients ψ σ 1 ...σ L ∈ C.
LetÔ be an operator acting on this tensor product Hilbert space H and letN = jn j be another operator with local operatorsn j : H σ −→ H σ fulfilling the commutation relations [n j ,n k ] = 0. We denote the ladder operators spanning the algebra of local operators byb ( †) j that obey canonical commutation relations b j ,b † k = δ j,k and = ± distinguishes between the commutator or anticommutator. Without loss of generality, we choose the spectrum of the local operatorsn j to be n j ∈ {0, 1, . . . , σ − 1} 2 . Let us assume furthermore thatÔ contains summands with ladder operatorsb Note that in this exampleN = jb † jbj is the operator counting the number of phonons and n f j measures the local fermion density. Next, we introduce a thermofield doubling of this Hilbert space. The new double Hilbert space H P B = H P ⊗ H B consists of two copies of the original Hilbert space, which we denote as the physical Hilbert space H P and the bath Hilbert space H B (see first arrow in Fig. 2). Correspondingly, we denote the density operatorsn P ;j andn B;j , which have exactly the same properties as the density operatorsn j in the original Hilbert space. In particular, the basis states |n P/B;1 ⊗ · · · ⊗ |n P/B;L ≡ |n P/B;1 · · · n P/B;L span a complete orthonormal basis of H P/B .
Here, we leave the framework of finite-temperature DMRG by considering the subspace P ⊂ H P B = H P ⊗ H B of the doubled system that is spanned by all states |n P ;1 , . . . , n P ;L ) = |n P ;1 , . . . , n P ;L P ⊗ |g(n P ;1 ), . . . , g(n P ;L ) B (9) = |n P ;1 , . . . , n P ;L , g(n P ;1 ), . . . , g(n P ;L ) P B , with n P ;j ∈ [0, σ − 1] and g(x) = σ − 1 − x (see second arrow in Fig. 2). Note that for convenience we have labeled the kets in the physical and bath system by subscripts and introduced rounded kets to indicate states in the subspace P ⊂ H P B , which depend only on a reduced number of coefficients n P ;1 , . . . , n P ;L . This subspace is contained in the subspace with N P + N B = (σ − 1) · L, i.e., (N P +N B ) |n P ;1 , . . . , n P ;L ) = (σ − 1) · L |n P ;1 , . . . , n P ;L ) , so that all states in the subspace P transform symmetrically under the action of the global U (1) symmetry generated byN P +N B . Furthermore, note that by counting the number of basis states spanning P it follows that dim H = dim P. Now, we define the map identifying states |ψ) ∈ P in the subspace of the doubled system with states |ψ ∈ H in the original Hilbert space as shown in Fig. 2. Since g(x) is invertible and dim P = dim H P = dim H, it follows that I is invertible. Next, we define the projected purified operatorÔ P P : P −→ P byÔ = IÔ P P I −1 .
AssumingÔ P P exists, this definition implies in particular that n 1 , . . . , n L |Ô|n 1 , . . . , n L = (n P ;1 , . . . , n P ;L |Ô P P |n P ;1 , . . . , n P ;L ) , that is, the matrix representations ofÔ andÔ P P in the local basis sets {|n 1 , . . . , n L } and {|n P ;1 , . . . , n P ;L )} are identical. We can, hence, work withÔ P P in the subspace P instead of O. In order to show thatÔ P P always exists, we construct it explicitly. For that purpose, we note that the above definition of P is equivalent to But this means that each operatorÔ P P has to satisfy Ô P P ,n P ;j +n B;j = 0 for all j ∈ {1, . . . , L} .
This motivates us to define balancing operatorsβ B;j /β † B;j : Since every operatorÔ P ⊗1 B acting non-trivially only on H P can be expressed as function of a product of ladder operatorsb P,j , we can thus map it to P through the transformationŝ and imposing the local gauge fixing conditions Eq. (15). By means of this transformation, which is shown graphically in Fig. 1, the local conservation laws Eq. (16) are fulfilled. There is also another way to introduce projected purified operators. We can define the projection operatorP = {nP;j} |n P ;1 , . . . , n P ;L ) (n P ;1 , . . . , n P ;L | and look for operators satisfyingPÔ P PP =Ô P P . Those operators are manifestly invariant under a projection into P and therefore, ignoring zero elements, have the same matrix elements in both H and P. Here the important observation is that restricting the ansatz class of states |ψ P B ∈ H P B to P, we have found a one-to-one mapping between H and P ⊂ H P B , and the states |ψ) =P |ψ P B transform under the global U (1) symmetry generated byN P +N B , obeying Eq. (15).
In the following, we explicitly derive the representation of states in P in terms of MPS and demonstrate the capability of the introduced U (1) symmetrization to improve the numerical efficiency of MPS calculations. For that purpose, we briefly recapitulate U (1)-invariant MPS before digging into the technical details of the projection. Consider a state |ψ as described in Eq. (7). Within the MPS formulation [31], the

U (1) Symmetries in Matrix-Product States
We refer to the matrix dimensions m j as bond dimensions. A compact representation of |ψ is then given by where neighboring matrices are contracted over their shared bond indices: . Commonly, these contractions are represented pictographically. Each tensor is drawn as a shape with as many legs attached to it as there are indices.
Then, contractions over shared indices are indicated by connecting the corresponding legs as shown in Fig. 3 for the case of a MPS. In order to exploit U (1) symmetries, let us consider a Hamilton operatorĤ : H −→ H of a system andN : H −→ H an operator generating a global U (1) symmetry ofĤ, i.e., with local density operatorsn j : H σ −→ H σ acting only on the jth lattice site. Since [Ĥ,N ] = 0, we can diagonalize both operatorsĤ andN in the same basis. Let this basis be spanned by {|N } withN |N = N |N as well as N |N = δ N,N . N is called the global quantum number of the state |N . A state |ψ ∈ H can now be expanded in terms of the simultaneous eigenstates |n 1 , . . . , n L ∈ H ofN with N = j n j and labels n j denoting the eigenvalues of the local operators 3n j : As a consequence of the Wigner-Eckart theorem, it can be shown [35,36] that the site tensors decompose according to with S . n j j; .
where we interpret in the following T . n j j; . π j−1 , Here, the indices α j are labeling irreducible representations of the U (1) symmetry on the bond spaces. Note that we introduce a dot over an index shifted towards the tensor . π j to indicate an ingoing index and a dot over an index shifted away from the tensor . π j to indicate an outgoing index. The orientation of the bond indices is also shown in Fig. 3. Hence, we can describe a state by its rank-5 site tensors M . n j j; . π j−1 ,   (15). Decomposition of the introduced auxiliary indexπ j−1 into irreducible representation of the local conservation law generated bŷ n P ;j +n B;j is sketched by double bonds γ j−1 → (ñ P ;j ,π j−1 ).

U (1)-Invariant Matrix-Product States with Bath Sites
The introduced mapping from an operator breaking a global U (1) symmetry to one conserving a U (1) symmetry (see Sec. 3) can be exploited to efficiently reduce the matrix sizes of MPS representations. The key observation is that, while purified states in the doubled Hilbert space in general have a huge redundancy that comes with additional gauge degrees of freedom, the projection into P fixes all these gauge degrees of freedom by the L local gauge constraints given in Eq. (15). Here, we discuss the implications on the projection of purified MPS into P and derive an important connection between the Schmidt decomposition of the purified states and the single-site reduced density matrix.
Let again |ψ ∈ H and consider its single-site representation with m j−1 |m j−1 = δ m j−1 ,m j−1 and m j |m j = δ m j ,m j . Following the previous considerations, we take this state representation into the subspace P of the enlarged Hilbert space H P B with N P + N B = (σ − 1) · L. We represent the MPS in H P B by interpreting the single-site representation of |ψ as a two-site representation in H P B , Then, we apply the projection into the subspace P by enforcing the local gauge condition Eq. (15). Pursuing those two steps at all sites j ∈ {1, . . . , L}, the resulting state representation is in the subspace P of the enlarged Hilbert space The resulting site tensors decompose under the global U (1) symmetry as A matrix factorization of the decomposed site tensors M . n P ;j , α j yields the MPS representation of |ψ) in the subspace P of the enlarged Hilbert space The MPS constructed in this way is shown in Fig. 4 and consists of alternating physical and bath sites, which are labeled by the physical and bath degrees of freedom n P ;j and n B;j , respectively. The delta function δ . n B;j ,g( . n P ;j ) in the last line of Eq. (32) is again the manifestation of the L gauge-fixing conditions imposed in Eq. (15). It motivates the introduction of the auxiliary bond labels η j enumerating the irreducible representations of each locally conserved quantity between the physical and bath sites. In this way the bond label γ j−1 can be decomposed into labels γ j−1 → (η j , a j−1 ) that need to fulfill η j + a j−1 = n P ;j + α j−1 . Note that we focus only on the labels for the symmetry blocks and -for convenience -in the following neglect the bond dimension m, which is part of the label π. From the local conservation laws and the gauge fixing, we can furthermore conclude that the bond label a j−1 has only one non-vanishing block with respect to the global U (1) symmetry, which is characterized by a quantum number (j − 1) · (σ − 1) ≡ α j−1 . Accordingly, there is only one non-vanishing block α j to the right of the bath site, which is characterized by a quantum number j · (σ − 1) ≡ α j . In tensor notation, this can be expressed by a reformulation of the local conservation laws at every site, introducing for brevity Therefore, we find that there is a unique decomposition of the auxiliary bond label γ j−1 = (η j , a j−1 ) given by identifying η j ≡ n P ;j and thus also a j−1 ≡ α j−1 . This can be summarized by decomposing the site tensors as γ j T . n P ;j j; .
which is exemplary shown in Fig. 5 and presumed from now on.

Connection to single-site reduced density matrix
The projected purification introduced above is closely related to the single-site reduced density matrix, which is the central object of the LBO method [13,40,42]. We consider the expectation value of the local density operators in the original Hilbert space written in terms of the singlesite reduced density matrixρ j = Tr k =jρ , n j = Tr j ρ jnj = n j n j |ρ jnj |n j = n j ρ n j ,n j n j .
Expanding the expectation value ofn P ;j in terms of the physical system's single-site reduced density matrixρ P ;j for states |ψ) ∈ P and a mixed-canonical MPS with center of orthogonality at the physical site j yields (n P,j ) = Tr P ;j ρ P ;jnP ;j = n P ;j (n P ;j |ρ P ;jnP ;j |n P ;j ) = n P ;j ρ n P ;j ,n P ;j n P ;j (36) = n P ;j ,n P ;j , n P ;j ,α j−1 , α j−1 n P ;j T . n P ;j j; . α j−1 ,( . n P ;j , . α j−1 ) δ n P ;j ,ñ P ;j * T . n P ;j j; . α j−1 ,( . n P ;j , . α j−1 ) δ n P ;j ,ñ P ;j = n P ;j n P ;j where we made use of the fact that the local symmetry generatorsn P ;j are one-dimensional representations of the local U (1) symmetry (see Fig. 6). From Eq. (14) it follows that Eq. (36) and Eq. (35) are completely equivalent so that ρ n j ,n j = ρ n P ;j ,n P ;j , and thus, comparing to Eq. (37), ρ n j ,n j = T . n P ;j j; .
We hence find that the single-site reduced density matrix of the physical part of P has the same diagonal elements as the original one. They are given by the trace over the absolute square of the symmetry blocks of the mixed-canonical site tensors. However, the symmetry n P ;j n P ;j n P ;j n P ;j · = T j T † j · δ n P ;j ,ñ P ;j δ n P ;j ,ñ P ;j δ α j−1 ,α j−1 α j−1 n P ;jαj−1ñP ;j n P ;j n P ;j Figure 6: Expectation value of the local density n P ;j , which by Eq. (39) can be directly related to the diagonal elements of the reduced single-site density matrix in the eigenbasiŝ n P ;j . conservation in P implies thatρ P ;j is diagonal whereasρ j in general is not. We can write the distance with respect to the 1-norm of these two operators by means of the mapping I: Here, we link to the LBO method, which expressesρ j in its eigenbasis (optimal modes) with diagonal elements w n j so that Let us now consider the Schmidt decomposition of a state |ψ at the auxiliary bond γ j−1 = (ñ P ;j ,α j−1 ). Because α j is fixed for every j, a block for a given n P ;j of a physical site can be decomposed individually to T . n P ;j j; .
The sum over the squared singular values is identified with the corresponding (diagonal) entry of the single-site reduced density matrix τ Λ j; . γ j−1 ,τ, . γ j−1 ,τ 2 = ρ n P ;j ,n P ;j .
Note that we implicitly accounted for all constraints arising from the projection into P and wrote the γ j−1 on the left only for completeness, as all α are fixed and the n P ;j is chosen. In Fig. 7, the argument is given diagrammatically. Truncating the singular values according to a certain threshold 0 < δ 1, so that n P ;j τ Λ n P ;j j;τ 2 < 1 − δ implies a rescaling of the diagonal elements of the single-site reduced density matrix (n P ;j |ρ P ;j |n P ;j ), which is governed by the decay of the singular values Λ n P ;j j,τ in each block. If we assume that the optimal modes ofρ j are truncated in the T j T † j α j−1 n P ;j γ j−1 n P ;j n P ;j ρ n P ;j ,n P ;j δ n P ;j ,ñ P ;j δ n P ;j ,ñ P ;j n P ;j Figure 7: For a given n P ;j , the diagonal entry of the single-site reduced density matrix is given by the singular values of the decomposed physical site. Note that we make extensive use of the tensor notation, in particular implicit deltas, which was introduced in [47].
same way, so that n j w n j < 1 − δ, we can compare this expression with Eq. (41). Then, using the invariance of the trace, a truncation of the bond index γ j−1 by means of the usual MPS truncation routine yields an equivalently precise approximation toρ j as the truncation occurring in the LBO. In addition, performing the truncation in the projected purified representation automatically favors those eigenvalues ofρ j that have the largest weight without the necessity of constructing the single-site reduced density matrix at all. This is an important improvement as it prevents the repeated constructions ofρ j in contrast to the LBO.

Characterization of numerical expenses
The previous considerations enable us to compare the numerical complexity of typical tensor contractions arising from the MPS representation of states |ψ) ∈ P with those of MPS representations without the expansion of the Hilbert space. At first, we point out again that due to the local conservation laws and the gauge fixing, the bond labels α j−1 , α j of the MPS site tensors T . n P ;j j; . α j−1 , . γ j−1 and T . n B;j j; . γ j−1 , . α j have only one non-vanishing entry; each of which is given by α j−1 = N j , α j = N j+1 with N j as defined above. Therefore, without truncation, the bond dimensions m j−1 , m j are identical to those of the site tensors M . n j j; . α j−1 , . α j representing the same state in the physical Hilbert space H only. There is no additional complexity arising from the representation of |ψ) ∈ P on these indices. Furthermore, without truncation the effective bond dimensions on the γ-bonds are given by m j;γ = n P ;j · min(m j−1 , m j ). In what follows, we analyze two truncation schemes on these bonds for states in the enlarged Hilbert space H P B . Thereafter, we discuss in which situations these yield a reduced numerical complexity of the most expensive operation during ground-state calculations, i.e., the application of a matrix-product operator (MPO) to a state.
A physically motivated truncation can be defined by exploiting Eq. (43) and discarding all single-site occupations ofρ j , whose sum is below a given threshold δ > 0. More precisely, let D ⊂ 0, · · · , n P ;j − 1 be a set for which n P ;j ∈D ρ n P ;j < 1 − δ. Sinceρ j is a reduced density matrix, its trace is normalized, and by sorting the diagonal elements such a set can always be defined. Then, all tensor blocks T . n P ;j j; . α j−1 , . α j−1 . n P ;j with n P ;j / ∈ D are discarded so that the total number of kept states on the auxiliary bond is bounded by m j;γ ≤ |D| min(m j−1 , m j ).
The physical interpretation is straightforward: All tensor blocks T . n P ;j that have a negligible single-site occupation T . n P ;j 2 =ρ n P ;j are discarded, i.e., empty modes do not contribute to the physics. However, we can give a tighter estimate by considering the explicit distribution of the singular values in each block.
Motivated by the numerical evidence that often the singular values decay exponentially in ground states of one-dimensional (1D) gaped systems [25,29,39], we assume such a decay in each block T . n P ;j j; . α j−1 , . α j−1 . n P ;j (n P ;j ∈ D). That means, in the decomposition shown in Fig. 7, Λ n P ;j j;τ = e −a n P ;j τ , for some a n P ;j > 0 and we abbreviated m j ≡ min(m j−1 , m j ). Note that n P ;j only specifies one block (due to the implicit δ n P ;j ,ñ P ;j ) and that we neglected the constant α j for brevity. Normalization to the single-site occupation yields Defining a n P ;j = − 1 2 log X n P ;j with 0 < X n P ;j < 1, we can rewrite Eq. (45) into Since δ ≤ ρ n P ;j ≤ 1 and m j ≥ 1, this equation has only one solution for X n P ;j in the given domain, even though there is no closed expression (see Fig. 8 for graphical solution at distinct pairs (ρ n P ;j , m j )). Therefore, we consider two limiting cases that yield upper and lower bounds on the decay of the singular values in each tensor block. The lower bound X n P ;j ,min is obtained through the intersection of the right-hand side with the horizontal axis and can be related to the limit m j 1: 0 = X n P ;j ,min (1 + ρ n P ;j ) − ρ n P ;j ⇒ X n P ;j ≥ X n P ;j ,min = ρ n P ;j An upper bound X n P ;j ,max can be established if the right-hand side of Eq. (46) is tangential to the left-hand side d dX n P ;j X m j +1 n P ;j X n P ;j ,max ! = 1 + ρ n P ;j ⇒ X n P ;j ≤ X n P ;j ,max = 1 + ρ n P ;j Combining both bounds, we find − 1 2m n P ;j log 1 + ρ n P ;j 1 + m j ≤ a n P ;j ≤ − 1 2 log ρ n P ;j which, by introducing normalization constants A n P ;j ,max/min , limits the decay of the singular values A n P ;j ,min X n P ;j ,min τ ≤ Λ n P ;j j;τ ≤ A n P ;j ,max X n P ;j ,max τ , and thus can be used to fix upper and lower bounds for the matrix dimensions required on the auxiliary bonds between physical and bath site. The normalization constants are determined from with η = min, max. We introduce a truncation threshold δ n P ;j for each block so that for singular values with τ ≤ m n P ;j ,η ≤ m j , we obtain ρ n P ;j − δ n P ;j ≥ A n P ;j ,η m n P ;j ,η τ =1 X n P ;j ,η τ = ρ n P ;j 1 − X n P ;j ,η m n P ;j ,η 1 − X n P ;j ,η m j ⇒ X n P ;j ,η m n P ;j ,η ≥ 1 − 1 − δ n P ;j ρ n P ;j 1 − X n P ;j ,η For this inequality to hold, we necessarily need ρ n P ;j − δ n P ;j ≥ 0, because A n P ;j ,η , X n P ;j ,η > 0. This is ensured by taking n P ;j ∈ D and choosing δ n P ;j = max( δ |D| , min n P ;j ∈D ρ n P ;j ) as truncation scheme. Then, taking the logarithm of both sides and solving for m n P ;j ,η , we divide by log X n P ;j ,η < 0 so that where we defined the truncation ratio R n P ;j = δ n P ;j ρ n P ;j ≤ 1. Imposing equality between the left and right side, we finally obtain an estimation for the upper and lower bounds of the required bond dimension m n P ;j ,η in each block. Introducing the relative change of the number of kept states F η (m j , ρ n P ;j ) = m n P ;j ,η m j , we show the bounds in Fig. 9 for varying m j and ρ n P ;j . For the upper bound there are two regimes: In the limit of small truncation ratio R n P ;j 1 we have F n P ;j ,max (m j , ρ n P ;j ) ≈ 1, whereas for R n P ;j → 1 there is a sharp drop towards zero. The transition regime between both asymptotics is governed by the physical bond dimension m j and shifts towards larger values of ρ n P ;j as m j increases. The lower bound exhibits a powerlaw decay over several magnitudes of ρ n P ;j and saturates towards one if m j is small (Fig. 9). Finally, from Fig. 8 we can deduce that if m j 1, the lower bound becomes an increasingly better approximation for the bond dimension m n P ;j ,j .
In summary, we found that for small physical bond dimension m j characterizing the approximation of the state without bath sites, the bond dimension m j,n P ;j between physical and auxiliary sites is of the order of |D |m j if m j is small (∼ O(1)) and D = n P ;j | ρ n P ;j > δ .
However, if m j 1, the relative value of the bond dimension m j,n P j per tensor block compared to m j mostly follows a power law in ρ n P ;j and quickly decays to zero. In this situation, the state can be efficiently approximated in the enlarged Hilbert space with a moderate growth of the bond dimension, given that the occupations of the single-site reduced density matrix ρ n P ;j decay fast enough.
In physical problems one is often faced with exponentially decaying occupations of ρ n P ;j [48,49]. Exemplary, we consider a typical, physical bond dimension m j = 100 and assume ρ n P ;j ∝ e −2n P ;j with a truncation threshold of δ = 10 −14 and take into consideration a local dimension of n P ;j = 21 (i.e., permit for 20 occupied states). We use the derived lower bound and obtain m j ≈ m j . This estimation relies on the assumption of strictly exponentially decaying singular values in each tensor block, which does not necessarily need to be the case in actual calculations. However, a relative growth in the overall bond dimension of O(1) was also found in our test calculations. Finally, we note that due to the rapid decrease of the lower bound derived above the total local dimension n P ;j is not a limiting factor in the first place as long as m j is large enough. In turn, the decay of the reduced single-site density matrix occupation strongly dictates the numerical expenses.
We close this section by demonstrating the numerical benefits of the above introduced enlargement of the Hilbert space and projection into the subspace P by considering the scaling of the most expensive calculation in a DMRG two-site ground-state search. This algorithm scales with the application of the MPO to the MPS and has dominating numerical expenses m 3 j · w j · n 2 P ;j if m j is sufficiently larger than w j . Assuming a typical growth factor 2 between the physical and bath sites, this operation is 8 times more expensive on these bonds than on the original bond between physical sites only. In order to benefit from the introduction of U (1)-invariant state representations in the first place, we therefore need to have a reasonably large local dimension n P ;j > √ 8, since for U (1)-invariant representations all local generators can be chosen as one-dimensional representations. Thus, n P ;j ≥ 3 already speeds up this contraction and the benefits will grow quadratically with larger n P ;j . We may also consider a decomposition of the MPO bond dimension w j due to the U (1) symmetry, which typically is of the order of 2 − 3 and thereby also generates an additional speed-up. Finally, we note that the system size is doubled, which could also be incorporated into the estimations. But this is only a constant factor of two and can be compensated easily by the quadratically growing expenses in the local dimension or the decomposition of the MPO bond dimension under the global symmetry.

Examples
In this section, we provide two simple examples, namely the Holstein model and the Hubbard model with superconducting (SC) terms. We concentrate on some technical details for the latter system, and present and discuss numerical results of the former.

Holstein Model
The Holstein model [50] is given bŷ in whichĉ ( †) j denotes spinless fermion annihilation (creation) operators,n f j =ĉ † jĉj the corresponding particle number operators, andb ( †) j the bosonic annihilation (creation) operators. The parameters of this model are the hopping amplitude t, the phonon frequency ω 0 , and the electron-phonon coupling γ. Here, the total number of spinless fermions jn f j is conserved, while the total number of phonons jb † jbj is not. Owing to the fermion-phonon interaction, the number of phonons per lattice site can become very large, rendering this model very challenging for DMRG, in particular in the charge-density wave (CDW) phase at half filling [6,7], for which we also present some numerical results.
We restore the conservation of the global phonon number by adding balancing operatorŝ β ( †) B;j , according to the procedure described in Sec. 3. The projected purified Hamilton operator then readŝ Note that the local phonon-density operators transform asb † jbj →b † P ;jbP ;j sinceβ † B;jβB;j = 1. Here, the last identity follows from the specific definition of the balancing operators in Eqs. (17) and (18).

Numerical results in the CDW phase
In order to illustrate the numerical properties of the mapping introduced in this paper, we performed calculations in the CDW phase of the half-filled Holstein model [7,40,51]. This phase is characterized by the formation of bound electron-phonon states (polarons) and a Fermi wave vector k F = π, i.e., in a physical image every second lattice site is occupied by a polaron. In the atomic limit t → 0, there is an analytic expression for the probability P ph (n j ) to measure n j phonons at occupied lattice sites j, which is given by Note that the excitation probabilities are given by the diagonal elements of the single-site reduced density matrix. Hence, they can be evaluated directly numerically. Another important quantity is the occupation w 0 of the optimal modes of the single-site reduced density matrix ρ j , which was already mentioned in Sec. 5.1. The optimal modes |d o are the eigenstates of ρ j and their occupations are the corresponding eigenvalueŝ As discussed elsewhere [13,40,42], these constitute an important measure for the quality of the approximation of the phonon states. In our framework, the full reduced single-site density matrix can be extracted directly from the projected purified state |ψ) in a mixed canonical representation when contracting physical and bath site tensors T n P/B;j over their auxiliary bond index γ j−1 (see Eq. (29): ρ j;n j ,n j = Tr k =j |ψ) (ψ| = Tr T n P ;j T n B;j † T n P ;j T n B;j , where we used the mapping I to identify n P ;j ≡ n j (see also Eq. (38)). For our calculations, we set ω 0/t = 1.0 and γ /t = 2.0 so that the model is in the CDW phase. In Fig. 10, the optimal modes of a system with L = 51 sites and N = 25 fermions are displayed for the ground-state and on an occupied lattice site (j = 25). The truncation was performed by allowing a maximum discarded weight of δ = 10 −14 per auxiliary bond while restricting the total bond dimension to m j ≤ 2000. The color-coded graphs correspond to calculations with different, maximally allowed total bond dimensions.
The immediate effect of the truncation on the auxiliary bonds between physical and bath site tensors can be seen as a suppression of the occupation w o of optimal modes when w o becomes small. Upon increasing the total bond dimension m j , the distribution w o (d o ) becomes stationary once m j > 1200. In the inset, the diagonal elements of the single-site reduced density matrix are shown as a function of m j and overlayed with the occupation probabilities P ph (n j ) (Eq. (56)) in the atomic limit. The discarded diagonal elements ofρ j can be deduced from the intersection of the vertical lines with the horizontal axis. Comparing the magnitude at which diagonal elements ofρ j are discarded as a function of m j to the plateaus of the optimal mode occupation in the main plot, we find a clear correspondence between both. This can be related to the discussion in Sec. 5.1, where we showed that w.r.t. to the 1-norm the quality of the approximation of the projected purified state is bounded by the occupation of the optimal modes ofρ j , which are not treated correctly. Thus, a scaling analysis in the bond dimension m j only is sufficient to obtain converged results for the phonon system. Finally, we find that, in accordance with the system being deep in the CDW phase, the diagonal elements ρ j;n j ,n j are already very close to the excitation probabilities P ph (n j ) in the atomic limit. Even though the bond dimensions m j ≤ 2000 may appear very large, the fact that we are able to exploit global U (1) symmetries for both the fermionic and bosonic system allows us to perform these calculations very efficiently. We also performed a finite-size scaling of the ground-state energy E δ as a function of the discarded weight per bond to prove the capability of our approach to deal with large system sizes. Here, we applied a scaling analysis in the numerical precision, tuning the maximal discarded weight per bond from δ = 10 −4 to δ = 10 −10 and extrapolated E 0 towards δ → 0. The number of lattice sites was increased from L = 51 sites up to L = 501 sites. In Fig. 11, we show the extrapolations and the scaling of the intensive energy density E 0/L as a function of 1 /L. We fit the ground-state energy densities as a function of the number of lattice sites using the ansatz Here, lim L→∞ E 0/L = ε ∞ is the extrapolated ground-state energy density in the thermodynamic limit yielding ε ∞ = −2.14628340 ± 4 · 10 −8 .
Note that the given uncertainty is obtained from propagating the errors of the scaling w.r.t. to the discarded weight per bond, which was done for each lattice size L. Since bond observables are evaluated with errors whose absolute values are bounded by the discarded weight per bond, this is a numerically exact error bound. Additionally, in the inset of Fig. 11, we plot the total CPU time of a ground-state search running until the convergence threshold Lδ with δ = 10 −8 for the relative change in the ground-state energy after a completed sweep was reached. Using two cores of an Intel R Xeon R Gold 6150 CPU @ 2.70GHz, the largest systems with L = 501 converged after ∼ 12 hours.
In this model, the pair creation contributions ∝ ∆ break the conservation of the global particle number conservation. We restore the corresponding global U (1) symmetry by adding balancing operatorsβ ( †) B;j,σ with σ =↑, ↓. The projected purified Hamiltonian now readŝ Again, the local density terms remain unchanged due toβ † B;j,σβB;j,σ = 1. Exploiting this representation, one of the authors studied the charge-degeneracy points of topologically superconducting islands coupled to normal leads [58][59][60][61].
In contrast to the Holstein model, here the projected purification acts on fermions. This causes a subtilty if the fermionic anticommutation relations are implemented in terms of Jordan-Wigner strings [62] as is usually done, either explicitly or implicitly. For instance, ifb ( †) j,↑ are annihilation (creation) operators of hardcore bosons at lattice j, then fermionic, bilinear operators can be written in terms of parity operatorsPb j,↑ aŝ The operator string k l=1Pb j+l,↑ is commonly referred to as Jordan-Wigner string and a consequence of the anticommutation relations. The problem here is that mapping such operator strings into the purified Hilbert space, one has to ensure that they act only in the physical Hilbert space. For instance, if the generation of the anticommutation relations is implemented in the MPS code itself, then typically such Jordan-Wigner strings are created automatically. If this is the case, their effect on the bath sites have to be canceled, which can be done by placing parity operators on bath sites inside the Jordan-Wigner string, for instance,

Conclusion
Numerically studying strongly correlated quantum many-body systems with a large number of local degrees of freedom is a challenging problem, in particular for tensor-network methods [6,7,[40][41][42]. In this paper we address the problem by introducing a mapping (projected purification) to construct artificial, global U (1) symmetries for models without a generic U (1) symmetry. For any given operator acting on a tensor-product Hilbert space H, we derived a construction scheme that generates its projected purified representation in a subspace of the thermofield doubling of H. We show that both operators can be identified with each other by an isomorphism, but the projected purified representation manifestly conserves global U (1) symmetries. Additionally, we derive a projected purified representation of MPS exploiting the fact that the isomorphism is obtained from a gauge fixing of the additional degrees of freedom introduced by the doubling. Here, the tensors representing the projected purified state can exploit the restored global U (1) symmetry which, for instance, immediately reduces the effective local dimension in each tensor block to 1 providing a significant speedup during numerical calculations when the local Hilbert space dimension is large. We characterize this representation and reveal an intimate relation between the Schmidt values of projected purified MPS and the single-site reduced density matrix that allows us to estimate the numerical expenses of our representation in comparison to calculations without symmetries. The mapping into a projected purified representation of operators and states is mostly independent of the underlying implementation. Thereby, it can be used without much effort with already existing toolkits, which we demonstrated by performing numerical calculations [63] on the one-dimensional Holstein model at half filling [7,10,51,64]. The large number of local degrees of freedom that have to be taken into account (we allow up to N ph = 63 phonons per lattice site) renders large scale calculations very challenging. We perform a finite-size scaling in the CDW phase taking into account a maximum number of L = 501 lattice sites while maintaining a high numerical precision and keeping up to m max = 2000 states per bond.
Due to the reduction of the effective local dimension of the MPS blocks, two-site solvers with a larger numerical complexity can be used [30,31,37,65], as we did in the ground-state calculations of the Holstein model. Therefore, the projected purification allows to apply twosite time-dependent variational principle (2TDVP) [32,66] as time evolution method to treat systems out of equilibrium. So far, existing methods to tackle such problems mostly [67] use time-evolving block decimation (TEBD) as time stepper, only, due to the high numerical costs when performing two-site updates on systems with a large number of local degrees of freedom [13]. However, TEBD typically requires a much smaller time step to achieve a certain precision, compared to 2TDVP [33]. We thus anticipate that using the presented mapping, out of equilibrium and finite-temperature calculations of such highly complicated systems can become cheaper, more reliable, and straight forward to realize. For instance, we expect this mapping to enable the efficient application of tensor-network algorithms to address questions about lattice electrons coupled to phonons out of equilibrium [22,68,69], numerically unbiased. Furthermore, our mapping is compatible with common MPO-based time-evolution methods, e.g., the aforementioned TEBD as well as the MPO W I,II methods [70]. Exhibiting a scaling of the numerical complexity that is at least quadratic in the local dimension [33], these time-evolution schemes should also benefit from taking operators to their projected purified representation.

A Object comparison between LBO and ppDMRG
Phys Bath π j−1 π j−1 π j = π j−1 + n P ;j + n B;j n P ;j + n B;j n P ;j n P ;j n B;j n j n P ;j n B;j n P ;j n P ;j + n B;j ω j−1 ω j = ω j−1 + (n P ;j − n P ;j ) + (n B;j − n B;j ) n P + n B n P ;j n P ;j n B;j n j n P ;j n B;j n P ;j n P ;j + n B;j π j−1 π j−1 π j = π j−1 + n P ;j + n B;j I : n B;j = g(n P ;j ) I † : n B;j = g(n P ;j ) I : n B;j = g(n P ;j ) I † : n B;j = g(n P ;j ) R † R R R † Id Id Figure 12: A tensor network representing a single site consisting of an MPS, an MPO, and the adjoint MPS. All tensors are split into several (virtual) objects in order to be rejoined to the tensors used in the LBO (ppDMRG) as highlighted by the red (blue) boxes that contain virtual objects. Note that equivalent bond labels do not indicate the same objects, but only an implicit δ between the, for brevity not shown, different indicies.
In this appendix, we aim to give an overview of the relationship between the objects used in the LBO and in the ppDMRG. Its main purpose is to support future discussions and developments. It is specifically not intended for implementation purposes, see Sec. 2. In Fig. 12, a complete sandwich MPS-MPO-MPS for a single site is shown. In order to show the connection between the LBO and the ppDMRG, all tensors are split into virtual objects that are subsequently rejoined in different fashions. On the one hand, the objects coming from the LBO (highlighted with red boxes) are mainly split vertically into parts "belonging" to the physical and the bath Hilbert space. On the other hand, the objects coming from the ppDMRG (highlighted with blue boxes) needed to be split horizontally so that they could be related to the different objects in the LBO. In particular, the identities containing the maps I and I † do not really appear within the ppDMRG.