On the Nature of Spinons and Holons - Elementary Excitations from Matrix Product States with Conserved Symmetries

We develop variational matrix product state (MPS) methods with symmetries to determine dispersion relations of one dimensional quantum lattices as a function of momentum and preset quantum number. We test our methods on the XXZ spin chain, the Hubbard model and a non-integrable extended Hubbard model, and determine the excitation spectra with a precision similar to the one of the ground state. The formulation in terms of quantum numbers makes the topological nature of spinons and holons very explicit. In addition, the method also enables an easy and efficient direct calculation of the necessary magnetic field or chemical potential required for a certain ground state magnetization or particle density.


I. INTRODUCTION
Matrix product state (MPS) 1-4 based methods such as DMRG, 5-7 TEBD 8 and VUMPS 9 have proven to be invaluable tools for simulating ground states of one dimensional quantum lattice models. By formulating those MPS methods in terms of manifolds and tangent spaces, 10 it has recently been shown that excitation spectra or dispersion relations as a function of momenta can readily be determined once the ground state is written in terms of a uniform (translation invariant) ground state. 11,12 Those tangent space methods extend the works ofÖstlund and Rommer,13 in which a slightly more limited ansatz was used. In this paper, we extend those tangent space methods to accommodate for U (1) symmetries, which are necessary to simulate quantum systems exhibiting a large amount of entanglement to good precision, such as e.g. the Hubbard model. The symmetric formulation further allows for targeting excitations with certain quantum numbers only, which greatly helps in disentangling rich excitation spectra of models which e.g. host several different types of elementary excitations.
In Sec. II we develop the theory of symmetric uniform MPS, while in Sec. III we introduce the necessary tools for formulating the excitation ansatz in the presence of symmetries, where we also generalize to multi site unit cells. This is done both for topologically trivial and nontrivial excitations which are domain wall like, such as spinons and holons. In Sec. IV, we demonstrate the usefulness of the methods by simulating excitation spectra of the integrable XXZ and Fermi Hubbard model, as well as the non-integrable extended Fermi Hubbard model. The use of symmetries makes the topological nontrivial nature of spinons and holons very clear and intuitive. Finally, we conclude with a summary and outlook in Sec. V.

II. SYMMETRIC UNIFORM MPS
We begin by defining properties for symmetric uniform MPS (suMPS), where we restrict the discussion to the case of abelian symmetries (e.g. Z n parity, or U (1) like particle number or magnetization). While symmetric tensor networks and the use of conserved quantities in tensor network algorithms have been addressed in numerous previous works 7, [14][15][16] , we reiterate here in detail their consistent use in the context of MPS in the thermodynamic limit.
In the following we closely use and follow notation and nomenclature of Ref. 9 and restate here only the most important concepts. For details we refer the reader to Ref. 9 (in particular Sec. II.A and II.E.). We consider a translation invariant uniform MPS in the thermodynamic limit in the mixed canonical representation Here A L , A R ∈ C D×d×D , with d the local Hilbert space dimension and D the MPS bond dimension, both describe the same state in different gauge representations and are related by the gauge transformation matrix C via where we have also defined the center site matrix A C . A L and A R then fulfill the left and right gauge constraints An N site unit cell uMPS, which is invariant under a translation over N sites, is described by N independent MPS matrices A(k), k = 1, . . . , N , which define the unit cell tensor where Σ n = (σ nN +1 , . . . , σ nN +N ) is a combined index. For an N site unit cell uMPS we then write where the integer n labels unit cells, not sites. Here we have not explicitly specified the gauge representations of the individual MPS matrices within the unit cell, but we will henceforth assume all states to be in the mixed canonical representation (1). Global symmetries of an MPS that is invariant under certain symmetry operations can be easily encoded in symmetry properties of the local MPS matrices A. For MPS with abelian symmetries this simply amounts to attaching quantum numbers to all indices appearing in tensor contractions and constraining the MPS matrices A to transform as irreducible representations of the global symmetry group. In practice, states on the physical and virtual level are therefore grouped into distinct quantum number sectors, and the matrices are of sparse block form. The action of the symmetry group then determines which combinations of quantum numbers are allowed, i.e. which blocks are non zero. We denote such symmetric MPS matrices as Here, a and b are quantum numbers labeling symmetry sectors on the virtual level and α, β label states within these sectors (degeneracy indices), while s denotes the quantum numbers associated with the local physical states σ, determined by a choice of numerical representation s(σ) (see below for examples), and ⊕ denotes the group action. 17 In the presence of several symmetries, the quantum numbers are multi valued. In the following we will often omit quantum number labels or degeneracy indices for better readability. To avoid ambiguity, we denote quantum numbers with Latin letters and physical/degeneracy indices with Greek letters. At times we also write -in a slight abuse of notation -A σ ab , where it is understood that b = a ⊕ s(σ).
In most cases the group action reduces to a simple (perhaps modular) addition/subtraction of properly defined quantum numbers, which we represent as (tuples of) rational numbers. In the following we will therefore denote a ⊕ b just as a + b, and the action with the inverse a ⊕b as a − b, where b ⊕b = 0. Furthermore, we regard s(σ) as a numerical representation of the physical state for quantum number arithmetic only, not necessarily as a strict group representation.
Without loss of generality, in (6) we have implicitly defined s and a as ingoing, and b as outgoing quantum numbers. 18 Upon concatenating symmetric MPS matrices, outgoing indices are then connected to ingoing indices only, and only sector blocks with matching quantum numbers are contracted. For example, a concatenation of two such matrices then yields Quantum states on finite systems of size L and with a certain quantum number Q can then be constructed by defining the quantum number of the left virtual boundary state to be zero, and the quantum number of the right virtual boundary state to be the desired quantum number Q. 19 Using symmetric MPS matrices (6), only basis states |σ 1 . . . σ L fulfilling the constraint j s(σ j ) = Q will then contribute to the overall state. This procedure is good practice and widely used in implementations of symmetry exploiting MPS algorithms. 4, 7 However, quantum numbers of ground states with U (1) symmetry in particular usually scale with the system size L (e.g. particle number N = L/2 at half filling), and the above scheme is therefore not scalable to the thermodynamic limit. In an equivalent, more homogeneous formulation, one can however demand both left and right virtual boundary vectors to have the same quantum number, 20 and to subtract the value of the desired quantum number densityq = Q/L as part of a modified group action for each individual MPS matrix Here we have defined a shifted quantum number for the physical index, which is the original quantum number s, offset by the desired overall densityq. This scheme is now easily scalable to the thermodynamic limit and we endow uniform MPS matrices with this modified group action (8) to obtain symmetric uniform MPS (suMPS) with well defined quantum number densities. The generalization to multi site unit cells is straightforward by endowing every matrix within the unit cell with the modified group action (8) and a consistent definition of the quantum number density. Consequently, the unit cell size N has to be chosen in accordance with the desired quantum number densityq. For example, a spin S = 1/2 suMPS with zero magnetization requires N even in order to host an equal number of up and down spins. Choosing a single site unit cell in this case results in a superposition of possible zero magnetization suMPS and hence a non injective MPS (i.e. the transfer matrix has multiple dominant eigenvalues with magnitude one).
We show concrete examples for shifted quantum numberss and required unit cell sizes N for spin S = 1/2 states with fixed magnetization densities m in Table I, and states of spinful electrons with fixed particle and magnetization densities (n, m) in Table II. 21 For an illustration of the conventional finite size scheme and the modified scheme for infinite systems, see Fig. 1. m − 0 1/6 −1/4 |↑ 1/2 1/2 1/3 3/4 |↓ −1/2 −1/2 −2/3 −1/4 N − 2 3 4 (1, 1/2) (0, 1/2) (1/3, 1/2) (−1/4, 5/8) It is worth noting that due to the quantum number density offset the shifted quantum numbers in general do not directly correspond to true quantum numbers of the symmetry group. This is because we have distributed the quantum number shift (to achieve a certain quantum number density) homogeneously over the unit cell, instead of applying a single offset at the unit cell boundary. Consequently, the shifted quantum numbers can be interpreted as a combination of the true quantum numbers of the symmetry group combined with quantum numbers of fractional applications of the unit cell translation operator (i.e. translations over n < N sites). The shifted quantum numbers therefore also bear information about the location within the unit cell. For example, in the 3 site unit cell suMPS with particle density n = 2/3 in Fig. 1, the sets of possible quantum numbers on the three bonds are mutually exclusive, i.e. on the first bond there are only integers, while on the second bond the shifted quantum numbers are integers shifted by +1/3, and finally on the third bond they are integers shifted by +2/3.

III. SYMMETRIC VARIATIONAL ANSATZ FOR ELEMENTARY EXCITATIONS
We now generalize the variational ansatz for low energy excitations presented Ref. 11 to multi site unit cells and abelian symmetries. As foundation we start from a variationally optimized suMPS ground state approximation |Ψ(A) of bond dimension D in mixed canonical representation, where for now we focus on single site unit cells and generalize to multi site unit cells later. A suitable suMPS ground state approximation can be obtained from e.g. a symmetric implementation of the algorithm presented in Ref. 9 using suMPS as described in Sec. II.
In a slight reformulation of the ansatz in Ref. 11 we write for a low-energy excitation with momentum p in the mixed canonical representation Here, the matrices left and right of B are in left and right canonical representation respectively. Furthermore, for a topologically trivial excitation A L andÃ R represent the same state, i.e. | Ψ(A)|Ψ(Ã) | = 1. For a topologically nontrivial excitation A L andÃ R represent different ground states (e.g. within the degenerate ground space in a symmetry broken phase), and the excitation is of domain wall type.
In the above mixed canonical representation, the parameterization of the perturbation matrix B in the left and right tangent space gauge 11 then reduces to Here V σ L and V σ R are the left and right null spaces of A σ L andÃ σ R respectively, i.e.
, and the matrices x L , x R ∈ C D×(d−1)D contain the (d − 1)D 2 variational degrees of freedom for the ansatz.
Without loss of generality we will henceforth use the left tangent space gauge (11a) and drop the subscripts L/R for V and x.
Variational approximations of excited states with a certain fixed momentum p can then be obtained from solving an effective eigenvalue problem of a (momentumdependent) effective Hamiltonian defined in the space of these variational parameters 11 where j = 1, . . . , (d − 1)D 2 and x denotes the vectorization of x. Notice that in contrast to the original formulation, in the above mixed canonical representation all necessary operations can now be performed without taking any (possibly) ill-conditioned inverses. For an N site unit cell ansatz, we again start from a variationally optimized (but here N site unit cell) ground state approximation |Ψ(A) in mixed canonical representation and introduce local perturbations in the form of local perturbation matrices B(k) on each site. We collect all of these single site contributions into one single unit cell perturbation matrix  Figure 1. Examples for possible quantum number sectors of MPS on a system of spinless fermions with fixed particle number. The lines represent possible paths, corresponding to valid sequences of available quantum number sectors. (a) Possible quantum number sectors in the conventional finite size representation of a 6 site system at half filling (i.e. particle density n = 1/2). Consequently the left and right boundary quantum numbers are 0 and 3 respectively, and the on site quantum numbers are s = 0, 1. (b) Possible quantum number sectors for an infinite system at 2/3 filling and a 3 site unit cell, with 3 quantum number sectors per bond. Here, the (shifted) on site quantum numbers ares = 1/3, −2/3. In typical low energy states, quantum number sectors corresponding to high fluctuations around the desired particle number density will in general be suppressed by small Schmidt values, which will be truncated away in a finite bond dimension MPS approximation, leaving a finite number of remaining sectors, even in infinite systems. Notice also the appearance of negative labels, as quantum numbers here are only defined up to an arbitrary global offset. In both cases we have marked in-and outgoing quantum numbers with arrows.
and write the full multi site unit cell ansatz with momentum 0 ≤ p < 2π N as Here again the integer n enumerates unit cells and we parameterize . . , N . The variational energy is then a quadratic function of the concatenation of all N parameter vectorizations x = k x(k), and multi site unit cell excited state approximations can be obtained from solving an effective eigenvalue problem of the same type as (12), but with a larger and more complex effective Hamiltonian H eff p (for an explicit construction of H eff p see Appendix A). Note that, contrary to blocking sites in e.g. a regular DMRG calculation, here the number of variational parameters scales linearly with N , enabling an efficient treatment of large unit cell sizes without sweeping.
By definition and construction of the variational ansatz and H eff p , the excitation energies e p obtained from (12) are (positive) energy differences to the extensive ground state energy E 0 . Hence, while E 0 is of O(L), the excitation energies e [j] p are of O(1). Likewise, while U (1) quantum numbers -like particle number or magnetizationare extensive for ground states, low lying excitations are characterized by quantum number differences of O(1) to the ground state. Popular examples for this are spin flip or few particle excitations.
In the context of the variational ansatz (10), such rela-tive differences of O(1) can be perfectly well understood as being caused by the single perturbation matrix B on top of the homogeneous, extensive ground state background generated by the MPS matrices A L andÃ R . We can therefore control and fix quantum numbers of excited states through the quantum numbers of B. More specifically, in the parameterization (11), V L and V R necessarily have the same symmetry sectors and quantum number labels as A L andÃ R . We can however attach a non-trivial quantum number q to the matrix x that contains the variational parameters, i.e. x is block (off-)diagonal and we write Note, that q is an (outgoing) quantum number associated with x itself, rather than a physical index. This is very intuitive, as V takes care of the homogeneous ground state density contribution on that site, while the quantum number q of x directly controls the quantum number difference of the excitation with respect to the ground state. For B we then have and we denote the generated symmetric excitation ansatz as |Φ p (B) . From the structure of (10) it is clear that the values of q have to be such that the outgoing quantum numbers of B match the ingoing quantum numbers ofÃ R . Depending on A L andÃ R , only certain values of q are therefore allowed.
For example, in a system of spinless fermions with s(σ) = 0, 1, a single particle excitation is characterized by q = 1, while a hole excitation corresponds to q = −1. Likewise, on a S = 1/2 spin chain with s(σ) = ±1/2, single magnons (spin flips) are characterized by q = ±1. This holds regardless of the ground state particle/magnetization density, which is encoded in V .
The above are typical examples for topologically trivial excitations, whereÃ R = A R , i.e. A L andÃ R have the same quantum number sectors, and the quantum number of the excitation is well defined. The quantum numbers of fractional excitations such as spinons and holons however necessarily require them to be of topologically non trivial nature (see below), whereÃ R = A R and the quantum numbers of A L andÃ R can in principle differ by an arbitrary offset. Just as the momentum p (cf. Ref. 11), the quantum number q of a topologically nontrivial excitation therefore seems to be completely arbitrary. Again, this is an artifact of open boundary conditions and fixing this ambiguity depends on the nature of the excitation (see Sec. IV A and IV C).

A. Spinons and Magnons in the S=1/2 XXZ Antiferromagnet
As a first prototypical example we study low energy excitations of the one dimensional S = 1/2 XXZ model in a magnetic field Here X, Y and Z are S = 1/2 spin operators, ∆ is the anisotropy parameter and h is an external magnetic field. The energies for the ground state and low energy excitations in the thermodynamic limit are known exactly. 22 We consider the antiferromagnetic case ∆ > 0. There the ground state in zero field has zero magnetization, and the elementary excitations are given by spinons. 22,23 In contrast to simple spin flip excitations (magnons) -which are integer spin excitations -spinons have fractional spin S = 1/2. 24 In the following we will demonstrate that in the context of the symmetric ansatz (14), spinon excitations must necessarily be of topologically nontrivial nature, i.e. they are domain wall like and cannot be created by a single (or few) local spin flips.
For a zero magnetization (m 0 = 0) suMPS ground state approximation, we require a unit cell size N = 2. We use a shifted quantum numbers(σ) = s(σ) = ±1/2, and the quantum number q of the excitation directly corresponds to the magnetization m of the excitation. 25 Without loss of generality we assume integer quantum numbers on even bonds, and half-integer quantum numbers on odd bonds. Consider the contribution For a topologically trivial excitation, the next unit cell starts again withÃ(1) R = A(1) R . The outgoing quantum numbers of V (2) and A(2) however are the same, which are in turn also equal to the ingoing quantum numbers of A(1), all of which are integers. This means that only integer values q are possible, such that both in-and outgoing quantum numbers of x [q] are integers. The local (single-mode approximation like) nature of the ansatz thus leads to the well known fact that excitations generated by localized spin flips can only generate integer spin excitations, i.e. magnons (where e.g. q = ±1 corresponds to single spin flips) The only possibility for q to be half-integer in order to generate a spinon, is for the unit cellÃ to the right of B to start with half-integer instead of integer quantum numbers. This can be achieved by using a translated unit cell forÃ, i.e.Ã(1) = A(2),Ã(2) = A(1), or |Ψ(Ã) = T |Ψ(A) with T the (single site) translation operator. Due to translation invariance of (17), |Ψ(Ã) is also a valid ground state approximation with the same energy as |Ψ(A) .
Indeed, in the gapped antiferromagnetic phase ∆ > 1 the exact ground state of (17) is twofold degenerate and spontaneously breaks translation invariance with finite staggered magnetization density m s = 1 L j (−1) j Z j = 0, and the above |Ψ(Ã) = |Ψ(A) are good ground state approximations. For −1 < ∆ ≤ 1 however the exact ground state is unique and the model is gapless. For finite D the above suMPS ground state approximations however artificially break translation invariance, such that still |Ψ(Ã) = |Ψ(A) are ground state approximations with the same variational energy and m s = 0 and we can use them to build spinon excitations. This symmetry is restored in the limit D → ∞, where m s → 0 and | Ψ(Ã)|Ψ(A) | → 1.
Colloquially speaking, the variational ansatz for a spinon is therefore obtained by "pulling" the 2 site unit cell ground state approximation apart by one site and inserting a variational perturbation matrix B carrying two half integer spins (one for the local physical S = 1/2 spin carried by V , and one for the spinon excitation carried by x) and the resulting ansatz is topologically nontrivial. Due to the 2 site unit cell, the momentum of the spinon is also restricted to 0 ≤ p ≤ π, i.e. to half of the first Brillouin zone, which is also consistent with Ref. 22 and 24. Here, generatingÃ from translating A also removes the above mentioned ambiguity of the excitation quantum number q, arising from the fact that A andÃ can in principle have arbitrary differences in the global quantum number offset. For a graphical representation of the construction of topologically trivial (magnon) and nontrivial (spinon) excitations, see    11 while an accurate representation of scattering states requires a more complicated ansatz, 26,27 whose symmetric formulation for multi site unit cells we leave for future work. Consequently, the elementary spinon branch is accurate to machine precision, while excitations in multi particle continua are only partially reproduced. Nevertheless, the low energy boundaries of these continua are still surprisingly well repro-duced, where accuracy however decreases quickly with higher particle number. An interesting advantage of the suMPS ansatz is that excitations very high up in the spectrum with high magnetizations (e.g. bound states) can now easily be targeted, whereas in the non-symmetric approach such states would be buried high up in a multitude of other states and targeting them would require more involved numerical procedures. This also enables a systematic estimation of unknown lower bounds of high up multi particle continua by targeting high magnetization excitations.

B. Magnetic Field for Ground State from Excitations in the XXZ Model
Apart from directly targeting and identifying excitations with certain quantum numbers, we can also use suMPS excitations to calculate the magnetic field h required for the ground state of (17) to have a certain magnetization density m 0 (or total magnetization M = Lm 0 ). This is possible, as the magnetic field term H M = j Z j commutes with the rest of the Hamiltonian H 0 , and changing h just changes the eigenenergies, but not the eigenstates of (17).
More  Table IV. Calculated magnetic fields h necessary for the ground state of (17) to have magnetization density m0 = 1/6. Shown are the values calculated from the data in Table III We demonstrate this method to calculate the magnetic field h required for a ground state with magnetization density m 0 = 1/6 in the antiferromagnetic cases ∆ = 0.5 and ∆ = 3. We start by obtaining the lowest energy state of H 0 with m 0 = 1/6, which requires a N = 3 site unit cell suMPS representation. Due to a finite bond dimension representation the obtained ground state approximation for ∆ = 0.5 also (artificially) breaks translation invariance, and for both values of ∆ there are three different and equally good ground state approximations which are related by single or two site translations. From these states we can construct several topologically nontrivial elementary excitations by combining the different translated unit cells left and right of B in an excitation ansatz with (fractional) quantum numbers q = ±1/3, ±2/3. There are 3 possibilities for each quantum number, totaling in 12 possible states. 28 See Appendix B for more details on obtaining excitations energies for well defined (fractional) magnetizations.
Without loss of generality we choose momentum p = 0 and calculate variational excitation energies e m 0,p with respect to H 0 with D = 200 and fractional magnetizations m = ±1/3, ±2/3 for ∆ = 0.5, 3. 29 The numerical results are given in Table III. From these values we further obtain the true excitation energies e m p with respect to H h from (21). These energies are known to be exactly zero 22 Fig. 4. This method is not restricted to determining necessary magnetic fields h for a certain ground state magnetization (density), but is generally applicable to all Hamiltonians which contain one (or more) generators of their global symmetries as a parameterized term. It is especially useful for models that are not exactly solvable, where otherwise one would have to perform a large number of ground state calculations in a grid search with small variations of h, while here h is calculated directly from a single ground state calculation followed by two (or few) excited state calculations.

C. Spinons, Holons and Electrons in the Fermi Hubbard Model
In the following two sections we consider as a second example the low energy spectrum of the (extended) Fermi Hubbard model where c σ,j , c † σ,j are creation and annihilation operators of electrons of spin σ on site j, n σ,j = c † σ,j c σ,j and n j = n ↑,j + n ↓,j are the particle number operators. Here, t is the hopping amplitude, U and V are the on site and nearest neighbor interactions respectively and µ is the chemical potential (we do not consider an external magnetic field).
Due to the phenomenon of spin charge separation, 31-33 the elementary excitations are fractionalized quasiparticles of either spin or charge alone, which cannot be constructed from the bare electrons, as those carry both spin and charge. Rather, electrons can in turn be interpreted as bound states of these elementary spin and charge excitations, known as spinons and holons respectively. Consequently, we use quantum number representation (n, m) of the local physical electronic states, with n the particle number and m the magnetization. Spinons and holons carry quantum numbers q s = (0, ±1/2) and q c = (±1, 0) and excitations with electronic quantum numbers q e = (±1, ±1/2) are therefore combinations of holons and spinons.
We first focus on the integrable case V = 0, where the ground state energy and elementary excitations are known exactly from Bethe ansatz, 30,34 at half filling (n 0 , m 0 ) = (1, 0), which requires a two site unit cell for a suMPS ground state approximation. Much similar to the case of the elementary spinon in the XXZ model in Sec. IV A, topologically trivial excitations can only carry electronic quantum numbers, which are unsuitable for creating elementary spinons or holons, and the elementary excitations are thus necessarily topologically nontrivial also in this case.
For the suMPS ground state approximation at half filling we use the corresponding shifted quantum numbers given in Table II. Even though the exact ground state is unique for all values of U > 0, the finite bond dimension suMPS ground state approximation again artificially breaks translation invariance, and the translated ground state |Ψ(Ã) = T |Ψ(A) = |Ψ(A) is an equally good ground state approximation. From these we can now construct topologically nontrivial excitations which allow for quantum numbers q c and q s and we can generate elementary spinons and holons that way. Conversely, topologically trivial excitations carry quantum numbers that are even combinations of q s and q c , for example spin-spin, charge-charge, or spin-charge excitations.
For the case N = 2 and half filling we calculate variational excitation energy dispersions for U = 5 and bond dimension D = 600 for several different quantum numbers. There, the elementary charge excitations are gapped, while elementary spin excitations are gapless. 30,34 The numerical results are shown in Fig. 5, together with exact results from Bethe ansatz, where we show pure spin excitations in the left plot and excitations with nonzero charge quantum numbers in the right plot. We find that the elementary spinon is reproduced up to an excellent accuracy of O(10 −6 ) by the lowest variational excitation branch with m = ±1/2. The higher up branches for integer and half integer magnetizations lie within the continuum of spin-spin scattering states which are not well captured by our ansatz. Nevertheless the variational energies reproduce the low end of this continuum surprisingly well.
Due to the elementary spinon being gapless, the exact elementary holon branch lies completely within the continua of scattering states of one charge and arbitrarily many spin excitations. Out of these, e.g. the chargespin-spin continuum (red area in Fig. 5) also contains excitations with quantum number q = (±1, 0) (charge + spin singlet or triplet with m z = 0) equal to q c . The variational excitation ansatz for this q tries to reproduce exactly these excitations, as due to the smaller bandwidth of the spinon branch they are at lower energies than the elementary charge excitations (shown as black solid lines in Fig. 5 on the right), except at p = 0, π, where the lower bound of the continuum is exactly the elementary charge branch. Around these momenta, the variational ansatz with q = (±1, 0) indeed yields the lowest energies and reproduces the elementary holon up to an accuracy of O(10 −5 ). Away from p = 0, π the same ansatz tries to reproduce a three particle scattering state, and energies for q = (±1, ±1/2) -which try to reproduce two particle p/π spin-charge scattering states -yield in fact slightly lower energies. A similar effect can be observed for q = (±2, 0) excitations, where there is a noticeable bend in the dispersion around p ≈ 0.15. There the ansatz tries to reproduce a spin-spin-charge-charge instead of a charge-charge scattering state. 35 Also, the exact elementary charge branch spans the entire Brillouin zone p ∈ [0, 2π), 30 while the momentum of a two-site ansatz is necessarily restricted to half of the Brillouin zone. However, it turns out that the two-site unit cell is required only by the spin quantum numbers, while half filling for charge quantum numbers alone could be achieved with a single site unit cell ansatz. Consequently, charge excitations are reproduced twice by this ansatz, with a relative shift of p = π over the entire Brillouin zone. For that reason we also draw the exact elementary charge branch from Bethe ansatz twice with the corresponding momentum shift in Fig. 5.
A possible way to remedy this fact is to realize that the variational ground state -even though not translation invariant under pure translations T -is invariant under a translation followed by a spin flip, i.e. under application of T F S , where T is the translation and F S is the spin flip operator. The entire Brillouin zone for charge excitations could therefore be recovered by using a restricted ansatz B(2) = e ip F S B(1), which indeed yields an eigenstate of T F S with eigenvalue e −ip where p ∈ [0, 2π) can be interpreted as quasi-momentum covering the full Brillouin zone.
Additional ways to only target the elementary charge branch within the charge-spin scattering continuum would be to either minimize the energy variance of the variational ansatz (see 27, especially Appendix 4), instead of the energy itself, or to distinguish excited states by higher conservation laws, 36,37 which could be introduced as artificial penalty terms into the Hamiltonian. The latter option however requires analytical knowledge about these higher conserved quantities and is unsuited for non-integrable models. We leave all these avenues to be explored in future studies.
Overall it is demonstrated that the new symmetric suMPS ansatz allows for an efficient separation of excitation sectors with different quantum numbers, which was possible in the original proposal in Ref. 11 only a posteriori and with great effort. For example, charge excitations only start appearing above the U dependent charge gap above a continuum of pure spin excitations. In the nonsymmetric original ansatz these excitations are next to impossible to single out or target, especially if the value of the charge gap is unknown.
In fact, even despite the elementary charge branch lying completely within multi-particle continua, the charge gap can still be calculated with the new symmetric ansatz to excellent precision of the same order as the ground state.

D. Spin to Charge Density Wave Phase Transition in the Half Filled Extended Fermi Hubbard Model
As a final example, we show the qualitative change of the low energy spectrum of the extended Fermi Hubbard model at half filling at the spin to charge density wave (SDW to CDW) transition for U, V > 0. 38,39 In particular we consider the case of U = 10 fixed and vary V around the critical point V c ≈ 5.13 of the first order SDW-CDW transition. While the charge excitations are always gapped in this parameter regime, the spin excita-  (23) in the non integrable case V = 0 at U = 10 on top of a half filled ground state (n0 = 1, m0 = 0). We show the first 10 lowest energies for each quantum number, represented by colored symbols. In addition we show spin-spin, charge-charge and spin-charge scattering continua constructed from the variational elementary spin and charge excitations as purple, green and red areas respectively. The top row shows spectra in the SDW phase V < Vc ≈ 5.13, while the mid and bottom rows show spectra for the CDW phase V > Vc. The SDW phase looks very similar to the integrable case in Fig. 5, while in the CDW case both spin and charge excitations are gapped and there is a multitude of additional isolated elementary (or bound state) branches. In the SDW phase -like in the integrable case -the lowest variational charge energies are only suboptimal approximations of multi-spin-charge scattering states; the charge-charge and spin-charge continua constructed from the variational energies are therefore not exact, but are kept for reference and marked with dashed boundaries. In the CDW phase the lowest excitation branches are isolated and the accuracy of variational energies -and consequently also of the multi particle continua -is expected to be excellent.
tions are gapless in the SDW phase (V < V c ) and become gapped in the CDW phase (V > V c ). We show results for variational excitation energy dispersions for V = {5, 5.1, 5.2, 5.5, 6, 6.5} and bond dimensions D = {400, 400, 200, 120, 100, 80} in Fig. 6 and plot the lowest 10 variational energies for various excitation quantum numbers. It can be seen that for V < V c in the SDW phase (top two panels in Fig. 5) the spin excitations are gapless and the dispersion looks very similar to the integrable case V = 0 in Sec. IV C. For V > V c in the CDW phase however, the spin excitations become gapped and the nature of the low energy spectrum changes completely: rather than one single elementary spin and one single elementary charge branch, there is now a multitude of isolated elementary (or bound state) excitation branches which lie below the multi particle scattering continua.
In particular, there are now two lowest excitation branches with spinon quantum numbers (n = 0, m = ±1/2) which have a level crossing at p = π/2. Likewise, there are also several isolated charge excitation branches with holon quantum numbers (n = ±1, m = 0). Here, the elementary holon is now also restricted to half of the first Brillouin zone, as the particle density n shows a strong dimerization in the ground state. While there are several level crossings between these charge branches for V slightly above V c , there remain two lowest crossing branches which become more and more separated from the rest with increasing V .
In addition, there are further isolated branches with spinon and holon quantum numbers, and also with electronic, or multi-particle spin or charge quantum numbers, which correspond to bound states. These branches lie completely or partially below the spin-spin, charge-charge and spin-charge multi-particle continua constructed from the lowest possible elementary spin or charge excitations (see colored areas in Fig. 6). For example, the lowest excitation branch with electronic quantum numbers (n = ±1, m = ±1/2) is on the same level as the elementary spin and charge branches and thus far below the spin-charge continuum, and even starts to lie completely below the elementary charge branch for V 6.

V. CONCLUSION
We present a formulation of the variational MPS ansatz for elementary excitations first proposed in Ref. 11 with conserved symmetries for multi site unit cells, where the computational cost and the number of variational parameters scales linearly in the number of sites N within the unit cell. The resulting ansatz allows for an efficient separation of the low energy excitation spectrum into certain desired quantum number sectors, which can be targeted individually. This is a great advantage over the original proposal, where an identification of different quantum numbers is only possible a posteriori, and there is no mechanism to target excitations with certain quantum numbers.
We show through the structure of the symmetric ansatz, that elementary excitations in the antiferromagnetic XXZ model (spinons) and in the Fermi Hubbard Model (spinons and holons) are necessarily of topologically nontrivial domain wall nature.
The performance of the proposed ansatz is demonstrated by calculating variational low energy dispersions with different quantum numbers for the antiferromagnetic XXZ model and the (extended) Fermi Hubbard model. In cases where exact Bethe ansatz solutions exist for comparison, the elementary spinon excitations are reproduced by the variational ansatz to excellent precision. In the gapped CDW phase of the (non-integrable) extended Fermi Hubbard model, we observe a large number of new bound states below the multi particle continua, which are not present in the gapless SDW phase. It would be interesting to explore the physical consequences of their appearance.
As the gapped elementary holon excitations in the (integrable) Fermi Hubbard model completely lie within a multi particle continuum with the same quantum numbers, the ansatz tries to reproduce lower lying states in this continuum instead, except around momentum p = 0, π, where the elementary excitation has the same energy as the lower boundary of the continuum. Possible ways to remedy this fact are discussed in Sec. IV C and are left to be explored in future studies. Despite this fact, the charge gap can however still be calculated with excellent precision of the same order as the ground state.
We further show that the symmetric ansatz can be used to calculate e.g. the magnetic field h required for a certain ground state magnetization m 0 of the antiferromagnetic XXZ model. As this strategy allows a direct calculation of h, involving a single variational ground state and two (or few) variational excitation calculations only, it is particularly useful for non-integrable models, where otherwise a large number of ground state calculations in a grid search with small variations of h are necessary. This procedure is applicable to all Hamiltonians, which contain generators of their global symmetries as parameterized terms.
The presented ansatz is also a perfect candidate for a more precise and efficient study of elementary excitations in two dimensional systems with topological order on cylinders, such as e.g. in Refs. 40  p onto some current vector x is necessary, which we describe in the following.
We divide the result x out = H eff p x in of one action of the effective Hamiltonian into individual contributions per site x(n) out , which are computed separately and combined at the end. Furthermore, we describe all terms in the space of B matrices, where B(n) σ in = V (n) σ L x(n) in and x(n) out = σ V (n) σ † L B(n) σ out . The individual contributions to H eff p can be derived by fixing the position of B(n) σ out and moving the positions of B(n) σ in and the two site Hamiltonian h n,n+1 (cf Ref. 11). Due to the gauge choice (11a) however, a good part of these terms are zero, namely those where B(n) σ in is strictly left of h n,n+1 and B(n) σ out , and those where B(n) σ out is strictly left of h n,n+1 and B(n) σ in . Below we give the remaining terms, collected and combined for efficient evaluation. Note that here the two site Hamiltonian has been offset by the energy density of the ground state (i.e. h → h − e 0 1 1 with e 0 = Ψ(A)|h|Ψ(A) ) in order to obtain positive energy differences to the ground state as eigenvalues of (12).
We follow the notation of Refs. 9 and 11 and write (x| and |x) for vectorizations of a D × D matrix x in the D 2 dimensional "double layer" virtual space, on which (mixed) transfer matrices, such as T A B = σĀ σ ⊗ B σ , or operator transfer matrices, such as O AB CD = σρµν O σρ µνĀ σBρ ⊗ C µ D ν (with O σρ µν = σρ|O|µν ), act. For better readability we also raise site indices to superscripts, e.g. A(n) σ L → A n,σ L and omit the tilde for A n,σ R . Furthermore we write T n L = T . In all expressions it is understood that N + 1 ≡ 1 and 0 ≡ N . We start by constructing quantities needed for terms, where B(n) in and B(n) out are on the same site. and and These quantities are independent of B in and p and can be precomputed as constants. They will also show up in other subsequent contributions. We now turn to quantities dependent on B in which have to be recalculated every time H eff p x in is invoked. We start with terms where B n in is right of B n out from within one unit cell and from all other unit cells Finally, we again collect all such contributions from B n in right of B n out into |B 1 R ) = e ipN |B 1 R ) |B n R ) = e ipN |B n R ) + |b n R ). (A9) These terms will be combined with (H n L | in the final contributions to B n out . Next we consider quantities where both h n,n+1 and B in are left of B out (A12) The inverses in (A8a) and (A11a) are to be understood as pseudo-inverses in the case p = 0 and q = 0 only (i.e. for excitations with zero momentum and the same quantum number(s) as the ground state) and can be fully inverted otherwise, as then the spectral radius of the transfer matrices is strictly smaller than 1.
We now have all the necessary quantities to compute B n where the Kronecker symbols δ n,N and δ n,1 ensure that the corresponding momentum factors only contribute in the cases n = N and n = 1 respectively. Here, the first line corresponds to contributions where B n in and B n out are on the same site, the second and third line where B n in is right of B n out , and the last line where B n in is left of B n out . For a graphical representation see Fig. 7.

Appendix B: Magnetization of Topologically Nontrivial Excitations from Translation Symmetry Broken Ground States in the XXZ Model
In this Appendix we elaborate on the topologically nontrivial elementary excitations with fractional magnetizations obtained from the threefold degenerate, translation symmetry broken ground states with magnetization density m 0 = 1/6 of (17) considered in Sec. IV B.
As a consequence of the broken translation symmetry, the single site ground state magnetizations m 0,n = 1/6, while N n (m 0,n − 1/6) = 0 still holds within one unit cell. Denote the suMPS unit cells for these three degenerate ground states as A 1 , A 2 and A 3 and assume |Ψ(A 3 ) = T |Ψ(A 2 ) = T 2 |Ψ(A 1 ) .
The resulting magnetizations of domain wall excitations created from combining these three different ground states are not exactly equal to the quantum numbers q. This is due to the fact, that the perturbation unit cell B consists of a superposition of different contributions from A i,L and A j,R (i = j), such that now n (m 0,n −m 0 ) = 0 within the perturbation unit cell. Consequently, B then has a magnetization different from N m 0 and the excitation carries a magnetization slightly perturbed away from q. These perturbations are usually of the form m 0,n −m 0 . For example, if we take the m 0,n to be the magnetizations of the A i,L unit cell, the q = 1/3 excitation then carries effective magnetization m = 1/3−(m 0,3 −1/6), while the q = −1/3 excitation carries m = −1/3 + (m 0,1 − 1/6).
However, the three possible excitations for each q (which are related by single or two site overall translations of the entire state) together have again a mean mag-netization of exactly m = q. For example, a q = 1/3 excitations can be generated from combining A 1,L with A 2,R , A 2,L with A 3,R or A 3,L with A 1,R . To remedy the above fact -which is once more an artifact of open boundary conditions and finite bond dimension -we therefore compute and average over all three excitation energies for each q. Exactly this has been done to obtain the values shown in Table III. Note that for topologically trivial excitations the magnetization is always well defined and precisely corresponds to m = q. This is because the same MPS ground state unit cell is used left and right of the perturbation matrix B.