Multiplons in the two-hole excitation spectra of the one-dimensional Hubbard model

Using the density-matrix renormalization group in combination with the Chebyshev polynomial expansion technique, we study the two-hole excitation spectrum of the one-dimensional Hubbard model in the entire filling range from the completely occupied band (n=2) down to half-filling (n=1). For strong interactions, the spectra reveal multiplon physics, i.e., relevant final states are characterized by two (doublon), three (triplon), four (quadruplon) and more holes, potentially forming stable compound objects or resonances with finite lifetime. These give rise to several satellites in the spectra with largely different spectral weights as well as to different continua. The complex multiplon phenomenology is analyzed by interpreting not only local and k-resolved two-hole spectra but also three- and four-hole spectra for the Hubbard model and by referring to effective low-energy models. In addition, a filter-operator technique is presented which allows to extract specific information on the final states at a given excitation energy. While multiplons composed of an odd number of holes do neither form stable compounds nor well-defined resonances unless a nearest-neighbor density interaction V is added to the Hamiltonian, the doublon and the quadruplon are well-defined resonances. The k-resolved four-hole spectrum at n=2 represents an interesting special case where a completely stable quadruplon turns into a resonance by merging with the doublon-doublon continuum at a critical wave vector. For all fillings with n>1, the doublon lifetime is strongly k-dependent and is even infinite at the Brillouin zone edges as demonstrated by k-resolved two-hole spectra. This can be traced back to the"hidden"charge-SU(2) symmetry of the model which is explicitly broken off half-filling and gives rise to a massive collective excitation, even for arbitrary higher-dimensional but bipartite lattices.


I. INTRODUCTION
The prime example for emergent phenomena in systems of elementary interacting particles is the formation of compound objects with new properties. Compound objects can be stable, i.e., described by exact eigenstates of the microscopic Hamiltonian. Common examples are the hydrogen atom, stable molecules, or mesons and baryons. In other cases, they are described by bound states interacting with a continuum. Those excitations acquire a finite lifetime and show up as resonances of finite width in various spectra, 1 such as Fano resonances 2 in the photo-absorption spectrum caused by autoionizing states of two-electron atoms, or the Auger effect in atoms and solids with inner-shell holes. 3,4 Recently, considerable attention has been devoted to repulsively bound pairs of particles 5 formed in a Bose-Hubbard model. 6,7 If the on-site Hubbard interaction U is sufficiently large compared to the width of the free Bloch band W , two bosons prepared on the same site move as a bound pair through the lattice. In an otherwise empty conduction band, this doublon is a completely stable compound object despite the fact that the Hubbard U is repulsive. This is a consequence of energy conservation which implies that there is no phase space available for a dissociation into two particles moving independently through the lattice. 8 The latter situation is described by states which belong to a scattering continuum at energies within twice the bandwidth 2W .
If there is, on the other hand, a finite concentration of bosons or fermions 9 in the conduction band, the doublon may decay by transferring its energy ∼ U in a high-order scattering process to several low-lying excitations of the continuum. This implies a finite lifetime of the doublon and poses an intricate many-body problem. The lifetime has been measured 10 for a well-controlled Fermi gas trapped in an optical lattice and was estimated theoretically for Fermi 10,11 and for Bose systems. 12 Studies of the real-time dynamics of an initially doubly occupied state |Ψ(t) = c † j↑ (t)c † j↓ (t)|0 , where |0 is the vacuum state or the few-particle ground state, have been done with exact-diagonalization methods. 13,14 For |0 being the ground state of the extended Hubbard model at half-filling, time-dependent density-matrix renormalization group (DMRG) 15,16 has been used to study the decay of a single doublon-holon pair 17,18 and of a large number of doublons. 19 In the context of condensed-matter systems, the idea of repulsively bound pairs goes back to the so-called Cini-Sawatzky theory developed in the late 1970s. 20,21 This addresses the energy-dependent cross section for corevalence-valence (CVV) Auger decay for a system with an initial core hole and a fully occupied conduction band described by a strongly correlated (Fermi-)Hubbard model H. Adopting a number of standard approximations (see Refs. 22,23 for a discussion), such as the sudden approximation and constant transition-matrix elements, and disregarding core-hole screening effects, the bare Auger spectrum I AES (ω) ∝ (−1/π) Im G 2−hole (ω + i0 + ) is given in terms of the local two-hole Green's function: Here, E 0 is the ground-state energy, ω is the two-hole excitation energy, i is the lattice site at which the two final-state holes are created, and |0 is the hole vacuum, i.e. the fully occupied band. Note that due to the core hole involved in the transition, the Auger process can be assumed as local to a good approximation. With a particle-hole transformation, the immediate relation to doublon dynamics becomes obvious: The Auger spectrum consists of a continuum of scattering states ("band-like part") related to the independent motion of the two final-state holes and, for sufficiently strong Coulomb repulsion U , of the so-called Cini-Sawatzky satellite at higher two-hole binding energies, which corresponds to repulsively bound two-hole states. It has a finite dispersion but zero width, so that the compound object of two holes moves through the lattice, but cannot decay. Calculations can be done analytically for the simple two-particle scattering problem. 20,21 However, beyond the case of a fully occupied conduction band, there is not much known about the two-hole excitation spectrum: The Cini-Sawatzky theory can be generated by a ladder approximation within diagrammatic perturbation theory. This idea has been used 24 to extend the theory to partially filled conduction bands, albeit in an essentially ad hoc way. Some progress was made by summing an extended diagram class which dominates in the low-hole-density limit 25,26 as well as by consideration of diagrams up second order in weak-coupling perturbation theory. 27 Except for straightforward exactdiagonalization studies, 28,29 which are plagued by strong finite-size artifacts, however, there has not been any systematic study of the two-hole excitation spectrum in the interesting regime of finite hole densities up to the halffilled conduction band where a strong Hubbard interaction is expected to be most effective.
With the present paper we discuss the two-hole spectrum in the entire filling range for the one-dimensional Hubbard model (including its extended variant) in the strong-coupling regime. To this end, we employ a matrixproduct state (MPS) implementation of DMRG 30 in combination with the Chebyshev polynomial expansion technique. [31][32][33] This method 34 provides us with essentially exact two-hole excitation spectra for lattices with several tens of sites. We also propose a filter-operator technique which can be implemented straightforwardly within the framework of the Chebyshev-DMRG technique and which is shown to be a very helpful tool to characterize the final states at a given excitation energy.
Our central goal is to study the fate of the concept of a repulsively bound pair of holes if it is allowed to interact with the states of the continuum of particle-hole excitations and how this manifests itself in the two-hole spectrum. This makes contact with several questions, for example: Is the doublon a well-defined resonance or even a stable compound object for a partially occupied band or does it decay immediately? How does its excitation energy, its dispersion and its weight change with the filling and can the excitation be traced continuously down to half-filling? Studying to those questions will put the clear physical picture of the Cini-Sawatzky model to a test.
In an effort to answer those questions, one soon recognizes that the two-hole spectrum is actually more complex than expected and exhibits several satellite features and continua. Those arise from heavier objects consisting of three and more holes. The question for the stability of these heavier multiplons arises immediately. To study the triplon and the quadruplon, in particular, it is helpful to consider higher-order spectroscopies as well, such as three-and four-hole spectroscopy for the Hubbard model or a two-hole spectroscopy for the hard-core boson model that captures the relevant low-energy physics close to complete filling.
The paper is organized as follows: The next section II briefly introduces the model and discusses Chebyshev-DMRG method. Results are presented and discussed in Sec. III. We start by recalling the physics of the completely filled band in Sec. III A. For the partially filled band, we first give an overview in Sec. III B before we discuss the main spectral features separately, namely the doublon (Sec. III C), the band-like part or two-hole continuum (Sec. III D), and the triplon and the quadruplon in Secs. III E and III F. In Sec. III G we introduce the new analysis tool based on filter operators and return to the question of the doublon lifetime in Sec. III H. Finally, a general and concluding discussion is given in Sec. IV.

II. MODEL AND METHOD
Using standard notations, the Hubbard model reads 35 We consider the model on a one-dimensional lattice of L sites with open boundary conditions and N electrons with L ≤ N ≤ 2L, i.e., for fillings n = N/L ranging from half-filling n = 1 to the limit of the completely occupied band n = 2. In Eq. (2), c † iσ creates an electron at lattice site i with spin projection σ =↑, ↓, and n iσ = c † iσ c iσ is the occupation-number operator. Summation over ordered pairs of nearest neighbors is denoted by ij . Furthermore, U ≥ 0 is the strength of the on-site Coulomb repulsion. Energy units are fixed by setting the nearestneighbor hopping t = 1, and the resulting width of the free tight-binding band is W = 4. It turns out as instructive to study the extended Hubbard model in addition. The Hamiltonian, includes an additional nearest-neighbor Coulomb interaction of strength V . DMRG has been previously employed to study its ground-state phase diagram. 36,37 One-particle spectral properties [38][39][40] as well as dynamical spin-and charge-correlation functions 41,42 have been computed using the dynamical DMRG method 36 and time-dependent DMRG. 15 Our goal is to compute the local two-particle spectral function which consists of the two-hole excitation spectrum with excitation energies ω ≤ 0, related to Auger-electron spectroscopy, and of the two-particle excitation spectrum with excitation energies ω ≥ 0, related to appearancepotential spectroscopy. 23 Here, |0, N is the ground state of H in the N -particle invariant subspace with the corresponding ground-state energy E (N ) 0 , while |n, N ± 2 denotes the n-th excited state with energy E (N ±2) n in the subspace two additional holes or particles, respectively. The corresponding chemical potentials µ (N ±1) are fixed by Particle-hole symmetry implies that the two-hole spectrum at filling n is related to the two-particle spectrum at filling 2 − n via It is therefore sufficient to study the two-hole and the two-particle spectrum for fillings at and above half-filling: 1 ≤ n ≤ 2. We employ DMRG in combination with the Chebyshev expansion technique. [31][32][33] The main idea is to expand the spectral function into Chebyshev polynomials T k (x) which are defined on the interval 0 ≤ x ≤ 1, either via a recursion relation, with T 0 (x) ≡ 1, T 1 (x) ≡ x, or with the weight function w (x) = 1/π √ 1 − x 2 from the orthogonality relation An application requires a proper rescaling of the spectrum of the Hamiltonian H to ensure that its rescaled eigenvalues have modulus less then unity. This can be most simply achieved with the linear transformation H = (H − b) /a where a = (E max − E min ) /2 and b = (E max + E min ) /2 and adding a small safety padding to avoid touching the edges. The extremal eigenvalues E min/max of H are obtained from a standard DMRG ground-state calculation in the N ± 2 subspace. A spectral function A(x) = 0 B 2 δ(x − H)B 1 0 of the rescaled frequency x, defined in terms of operators B 1 and B 2 , is obtained as Using Eq. (10), the moments µ k can be expressed as and are calculated as µ k = 0|B 2 |t k from states |t k which, by making use of the recursion relation Eq. (9), can be constructed recursively, starting from t 0 = B 1 0 and t 1 = H t 0 . A final back-scaling yields the desired A(ω). Naive truncation of the series in Eq. (11) at k = M < ∞ causes so-called Gibbs oscillations close to sharp features in the spectral function. 31 Those can be attenuated by several orders of magnitude when multiplying the moments µ k by correction factors g . This amounts to a convolution of A(x) with a kernel and results in the kernel-polynomial approximation The two main choices for g (M ) k are the Lorentz and the Jackson kernel 31 corresponding to Lorentzian or Gaussian broadening of the δ-peaks, respectively. While the Lorentz kernel preserves the analytical properties of the spectral function, the Jackson kernel offers a somewhat better resolution. When working with matrix-product states, the effort in calculating the moments becomes quite significant, and we therefore opt for the latter throughout this work in order to obtain more spectral information for a given value of M . We have not seen any need to apply linear prediction 30,43,44 in the present study.
We apply the Chebyshev expansion method to exact many-body states for moderately large Hilbert spaces and, in case of large systems, combine it with the concept of matrix-product states 30 as suggested in Ref. 34. An MPS is obtained by a matrix-product ansatz for the coefficients of the many-body wave function in the product basis of states { s i } spanning the local Hilbert space at a site i: For the Hubbard model, s i = 0 , ↑ , ↓ ↑↓ . This ansatz allows an efficient truncation of the wave function by truncating the individual matrices. Furthermore, many operations can be applied locally with subsequent sweeps through the chain until convergence is reached. However, the MPS representation also comes at a cost, as the application of H and the sum of two states in Eq. (13) can usually be done in an approximate way only. Namely, to find an optimal representation of the left-hand side of Eq. (13), we perform a variational compression, 30 sweeping along the chain and optimizing one site at a time until t k − (2 H t k−1 − t k−2 ) 2 < with a given tolerance . Throughout this work = 10 −4 is used. Note that it is possible to set up two environments for the overlaps t k 2 H t k−1 and t k t k−2 , and impose the above condition locally without the need of carrying out the subtractions one at a time. As an initial guess to start the sweep we use the Chebyshev vector from the previous iteration. If convergence is not reached after four half-sweeps, we switch to a two-site algorithm which allows both to escape a local minimum and to dynamically increase the matrix dimensions, while keeping track of conserved quantum numbers.
The application of a local operator B 1 = c i↑ c i↓ to the ground state 0 creates a local perturbation which propagates through the lattice when repeatedly applying the Hamiltonian. The entanglement entropy therefore grows with the iteration index k, i.e., the state t k in general requires a higher bond dimension than t k−1 for a correct representation of the wave function. We typically observe a logarithmic increase of the entanglement entropy with k, consistent with the findings of Ref. 45. The maximal iteration depth is hence limited in DMRG as, assuming a fixed error, one eventually either runs out of memory or the procedure slows down beyond practicability.
The Chebyshev method has an approximately uniform energy resolution δE given by δE ∼ ∆E/M where ∆E = E max − E min is the total bandwidth. 31 In order to be able to compare the spectra for various fillings and pa-rameters, we therefore fix δE as a measure of resolution rather than M in all calculations.

A. Completely filled band
Let us start the discussion of the results with the simple case of the completely filled band: n = 2. Creating two holes with opposite spins constitutes a two-body problem and a closed analytical expression for the excitation spectrum is readily derived. 20,21,24 Alternatively, one may apply the Chebyshev expansion method using exact states for rather large systems, since the Hilbertspace dimension is merely L 2 . Some numerical results for L = 1000 and periodic boundary conditions are shown in Fig. 1.
The two-hole spectra for different U can be understood as follows: If U = 0, one can straightforwardly evaluate Eq. (5) by Fourier transformation to reciprocal space. This yields where ρ 0 (x) = L −1 k δ(x − ε(k)) and ε(k) = −2t cos(k), i.e., the two-hole spectrum is given by the self-convolution of the non-interacting density of states. 46 The self-convolution of the one-dimensional density of states ρ 0 (x) is just given by the two-dimensional one with bandwidth 2W = 8t (see black line). The two conduction-band holes move independently of each other through the lattice.
Upon increasing U , a sharp satellite emerges at the lower edge of this "band-like part" of width 2W and moves to higher binding energies. 20,21 For U t it is approximately separated by U from the barycenter of the band-like part located at ω = −4t. This corresponds to the formation of a stable two-hole bound state which we will call "doublon" in the rest of the paper. Since the one-hole spectrum has a finite support, a decay of the doublon into independent holes is prohibited by energy conservation for large U . One dimension is in fact special in that there is a bound two-hole state outside the continuum for any U > 0, so that momentum conservation must be considered as well to understand its stability. 8 This is very similar to the Bose case discussed in Ref. 5.
The finite width of the satellite reflects the fact that the doublon is itinerant. For strong U the propagation of a single doublon is perturbatively described by a secondorder hopping process and is captured by the effective Hamiltonian (see also Eq. (27) below) where d † i = c † i↑ c † i↓ is a bosonic doublon creator and n d i = d † i d i . Hence, in the strong-coupling limit the satellite has the shape of the non-interacting density of states with a rescaled bandwidth of 8t 2 /U . With increasing U , the doublon satellite gains more and more spectral weight at the expense of the band-like part and eventually results in a δ-peak for U = ∞, corresponding to a completely immobile doublon.
In the strong-U limit, the physics can be visualized to some extent by typical electronic configurations contributing to the different spectral features. Fig. 2 gives a sketch of the different configurations contributing to the doublon satellite "D" and to the band-like part "B".

B. Partially filled band: Overview
We now turn to the case of a conduction band with a finite hole density: n < 2. This changes the situation completely as the two additional holes placed with opposite spins at site i interact with a thermodynamically relevant fraction of holes in the incompletely filled band. From the technical point of view, the computation of the spectrum transmutes into a severe many-body problem: A comprehensive and reliable study, however, is possible when tolerating a non-zero low-energy cutoff in the form of a finite spectral resolution δE. Using the Chebyshev technique combined with matrix-product states, we have been able to study the excitation spectra in the entire filling range 1 ≤ n ≤ 2 for the Hubbard model on lattices as large as L = 60 at a constant energy resolution δE ≈ 0.2t. With δE 0.2t one would start to resolve finite-size artifacts. Calculations have been performed for systems with open boundary conditions and choosing i as a central site in the chain [see Eqs. (5) and (6)]. Figure 3 provides an overview and shows the filling dependence of the spectrum at U = 6t, chosen such that for n = 2, the doublon satellite is clearly separated from the band-like part. Several distinct effects arise at finite hole density, which we discuss now in turn: First of all, the overall intensity of the two-hole spectrum in the range ω < 0 is seen to decrease with decreasing filling. This is a consequence of the sum rule for the double occupancy, which is straightforwardly derived from Eq. (5), and of the fact that n i↑ n i↓ is decreasing with n. The doublon satellite ("D") is still present with strong weight at U = 6t but its support broadens significantly with decreasing n. The satellite can be traced down to fillings close to half-filling. Right at n = 1 it vanishes completely. An extended discussion is given in Sec. III C.
A continuation of the "band-like part" ("B") of the two-hole spectrum is also clearly visible, at least down to n ≈ 1.5. When the filling approaches half-filling, however, the weight of the band-like part becomes very small, and the structure is hardly visible at low energies, say −2 ω < 0, for fillings 1 < n < 1.2. Note, however, that there is finite spectral weight for ω → 0 at any filling n > 1; only at half-filling, n = 1, is the spectrum gapped.
As soon as n < 2, the spectrum exhibits an additional structure, marked as "T" in Fig. 3. For fillings n → 2 it shows up as a tiny peak at high binding energies. It gains some spectral weight upon decreasing n and moves towards ω = 0, until it completely merges with the doublon satellite. At the same time, another peak can be seen emerging within the spectral range of the doublon satellite (seen best for n = 1.6 at ω ≈ −8.7t). It moves towards ω = 0 at the same rate and also belongs to "T". The structure "T" will be interpreted as a "triplon" in Sec. III E. The physics close to and at half-filling is complicated by the appearance of yet another structure "Q" which emerges for n 1.5 at the highest binding energies and gains more and more weight as n approaches half-filling. This feature will be addressed in Sec. III F and interpreted as a "quadruplon".
At half-filling the system is a Mott insulator. The gap in the one-particle spectral function is accessible to the Bethe-ansatz method 47 and can be computed from the exact expression 48 At U = 6t it amounts to ∆ = µ + − µ − ≈ 2.89t. As Fig. 3 demonstrates, the two-particle spectral function is gapped as well, but by twice the one-particle gap 2∆ ≈ 5.86t (the numerical value is obtained from the groundstate energies, and the slight discrepancy compared to twice of the Bethe-ansatz value is due to the finite system size). Note that the chemical potential jumps from µ = µ − for n 1 to µ = µ + for n 1 and is not uniquely defined at half-filling. 49 The two-hole and the two-particle spectrum for n = 1 in Fig. 3 have been calculated for the choice µ = (µ + + µ − )/2 = U/2 = 3t, which implies a symmetric spectrum A 2 (−ω) = A 2 (ω). If µ = µ + had been chosen, the n = 1 spectrum would have evolved continuously from the spectrum for n > 1. Hence, the single broad peak in the two-hole spectrum at n = 1, which takes the whole spectral weight, smoothly connects to the quadruplon peak "Q" at higher fillings.
For fillings below half-filling, n < 1, the two-hole spectrum is related to the two-particle spectrum for n > 1 via the relation Eq. (8). As is seen for ω > 0 in Fig. 3, it consists of a single peak, smoothly connected to "Q" which just continues to lose spectral weight until it vanishes for n → 0. The vanishing of the weight is related to sum rule for the two-particle spectrum (20) which is straightforwardly derived from Eq. (6), and to the fact that the average number of empty sites goes to zero for n → 2, i.e., the phase space for creating two additional particles at a site i vanishes.
Let us also mention that the difference between the total weights of the two-hole and the two-particle spectrum is given by the filling: This is a direct consequence of Eqs. (18) and (20).

C. The doublon
For fillings close to the simple n = 2 limit, the physical interpretation of the doublon peak "D" as a repulsively bound pair of holes is still self-evident. With decreasing n, and for fillings close to half-filling in particular, the validity of this picture must be questioned seriously.
First of all, as soon as n < 2, the doublon starts to interact with a continuum of electron-hole excitations, such that the corresponding peak should actually be seen as a resonance with a finite lifetime rather than an exact eigenstate of the Hamiltonian (for a given wave vector). This will be discussed in detail in Sec. III H.
In fact, the doublon peak broadens with decreasing n. Its support widens from a value slightly larger than ∆ω = 8t 2 /U = 4t/3, which is the perturbative result for U → ∞ at n = 2, to ∆ω 4t for fillings close to, but above half-filling.
To understand the main reason for this broadening, consider the spectrum for a system with N = L + 1 particles, i.e., for a filling slightly above half-filling n = 1. Fig. 4 provides a sketch. The support of the doublon peak is given by the spectral range of the final states. As indicated in the sketch, the dominating final-state configurations are characterized by a hole moving through a lattice of half-filled sites. This is a standard manybody problem, intensively discussed, e.g., in the context of the Mott-Hubbard metal-insulator transition: 50 Double occupancies are effectively suppressed at strong U , well-formed local magnetic moments emerge and develop strong antiferromagnetic correlations due to the Anderson super-exchange mechanism. While the effective coupling constant J = 4t 2 /U of the latter is derived in second-order-in-t perturbation theory and is small, the motion of the hole takes places on the large energy scale t. Neglecting J altogether, already provides the right width of about 4t for the doublon peak, consistent with the corresponding spectrum displayed in Fig. 3. Let us mention that the arguments can be formalized by considering the standard mapping of the low-energy sector of the Hubbard model for strong U onto the t-J model. 51 Close to half-filling it reads: where S i = 1 2 σσ c † iσ σ σσ c iσ is the spin at site i, σ the vector of Pauli matrices and where c iσ is the annihilator in the subspace with no double occupancies. Neglecting J yields the t model which exactly maps 52 onto the J KL = ∞ Kondo lattice. For a single hole in the t model or, equivalently, a single electron in the Kondo lattice, the ground state is a ferromagnet as is known from Nagaoka's work 53 or, equivalently, can be proven 54 for the Kondo lattice. The width 4t is obvious then. For the present case with finite U , a somewhat larger width is expected.
Another helpful observation is that the final states of the two-hole spectroscopy (see Fig. 2b) are the same as the final states of the one-hole spectroscopy, i.e., photoemission, but from an exactly half-filled conduction band with N = L electrons. The (raw) photoemission spectrum is given by the one-hole spectral function . (23) This shows that the support of the doublon peak must be the same as the support of the photoemission spectrum. The latter has indeed a width of ∆ω 4t at U = 6t, as can be read off from the DMRG results shown in Fig.  19 of Appendix A. As the matrix elements for the twohole spectroscopy are different, however, the shape of the doublon satellite is different and, at least on this level, it is difficult to make contact with the well-known spinon and holon excitations that govern the one-hole spectrum.
We conclude that the concept of a doublon as a repulsively bound pair of holes breaks down for fillings close to half-filling. The main difference is that the doublon moves in a background of singly occupied sites and explores antiferromagnetic correlations. This strongly affects its mobility. For n → 1 the doublon propagates on the energy scale of about 4t, opposed to the energy scale 8t 2 /U set by double-hopping processes in the highdensity limit. The strong decrease of the weight of the doublon peak with decreasing n is related to the decreasing density of doubly occupied sites in the respective initial ground state. The shift of the position of the doublon satellite is closely related to the filling dependence of the band-like part.

D. The band-like part
The band-like part "B" of the two-hole spectrum results from the continuum of unbound two-hole final states. For n = 2 and U = 0, this is given by the self-convolution of the density of states ρ 0 (x). With increasing U it quickly loses weight in favor of the doublon satellite (see Fig. 2) and becomes an almost featureless structure for strong U whose shape still loosely resembles the self-convolution of the density of states, but without the strong van Hove singularities. Its weight is already low for U = 6t > W , see Fig. 1.
For fillings close to half-filling, this scenario still holds to a good approximation, but with the density of states replaced by its occupied part only (for a given filling). Since the occupied fraction shrinks with decreasing n, the weight of the band-like part must diminish and its spectral support narrow down, with the barycenter shifting upwards. This effect is clearly visible in Fig. 3 for high fillings.
The barycenter of the band-like part also fixes the position of the doublon satellite which is found at an energy of about U below. Assuming a rectangular density of states with bandwidth W = 4t, we have ω 0 = −(n/2)W for the barycenter and ∆ω = nW for the width of the band-like part. This simple model roughly explains the trend seen for fillings n 1.6. For example, the model predicts ω 0 = −3.2t for n = 1.6 resulting in ω = −(n/2)W − U = −9.2 for the satellite position. This should be compared with the barycenter of the satellite at ω ≈ −8.75t in Fig. 3 (obtained by integrating the data).
For fillings n 1.5, however, ω 0 shifts much faster with decreasing n and the model breaks down (note that there is a finite band-like part for any n > 1). In this filling range, one may consider the two additional holes in the final state as moving independently of each other, but strongly interacting with the high hole density in the ground state.
This can be very roughly described in the following way: We replace the non-interaction density of states ρ 0 by an interacting density of states ρ which is assumed to consist of two Hubbard peaks, separated in energy by U . For n ≥ 1 and following the standard theory 55 of filling-dependent spectral-weight transfer in the strongcoupling limit, the lower Hubbard band has the weight 2−n and is fully occupied, while the upper Hubbard band is partially occupied with the occupied fraction having the weight 2n − 2 (and the unoccupied fraction having weight 2 − n). In the simple model the band-like part is then proportional to the self-convolution of the lowenergy part of the occupied density of states ρ, i.e., of the occupied part of the upper Hubbard band. Assuming that the shape of the latter is n-independent, we find This roughly explains the filling dependence of the doublon satellite for n 1.6t. For example, we have ω 0 = −0.8t for n = 1.2, resulting in ω = −(n−1)W −U = −6.8t for the satellite position. This compares well with the DMRG results in Fig. 3 where we find ω ≈ −6.3t. For n → 1, in particular, the model predicts a vanishing band-like part with ω 0 → 0, and the resulting position of the doublon satellite at ω ≈ −U agrees well with the data.
Although the model is surely oversimplified, it illustrates that the filling-dependent spectral-weight transfer is important for the location of the continuum of unbound two-hole final states, and thus for the doublon satellite as well.

E. The triplon
The most apparent deviations from the standard Cini-Sawatzky theory consist in the appearance of the new structures "T" and "Q". Let us first concentrate on "T" which already shows up for low hole-densities n 2, a limit which has been studied to some extent using the ladder approximation in the hole-hole channel. 25,26 When acting with c i↑ c i↓ on the ground state, one may distinguish between contributions to the spectrum resulting from two different types of ground-state configurations: In the case of an excitation that takes place in a hole-free region, one expects essentially the same shape of the spectrum as in the n = 2 case. For n → 2 this is the dominant contribution. If the excitation takes place at a site i next to an already existing hole, however, a final state with three holes on nearest-neighboring sites is created. This case is schematically illustrated with Fig.  5.
Here, the question arises whether the three holes can form a new bound object, a "triplon". Since this question addresses the excited states one projects onto in Eq. (5), rather than the ground state, we first study the threehole excitation spectrum for the completely filled band (n = 2). This involves the same type of final states while the ground state is trivial, and we are left with a much simpler three-body problem only. Since the Hilbert-space dimension is approximately L 3 /2, we can use the Chebyshev expansion method with exact states.
The result for L = 500 sites is shown in Fig. 6. One observes a similar general pattern as in the two-hole case:  For U = 0 the spectrum consists of a three-fold selfconvolution of ρ 0 (ω) with a spectral width given by 3W = 12t, while for finite U > 0, there is a band-like part and a split-off correlation satellite.
For finite but strong U , the satellite in A 3−hole,σ (ω) has a slightly asymmetric weight distribution due to the presence of the band-like part. Apart from that, its shape is given by a convolution of ρ 0 (ω) (the single hole, width W ) with the same density of states, but rescaled (the doublon, width ∆ω = 8t 2 /U ), see the inset in Fig. 6. This clearly indicates that the doublon and the hole propagate independently and do not form a bound state.
This should also hold for the two-hole spectrum (Fig.  3): The structure "T" does not represent a stable triplon, but is rather a doublon-hole continuum superimposed on the doublon satellite. The filter-operator analysis in Sec. III G and further arguments given below indeed support this interpretation.
Compared to the three-hole spectrum, "T" in Fig. 3 has smaller spectral weight, but the same shape. We would like to stress that the two peaks belonging to "T" are simply a consequence of the van Hove singularities in one dimension, separated by about 16t 2 /U ≈ 2.67t for U = 6t, and should not be mistaken for separate resonances. Hence, they will not be observed for higherdimensional lattices.
The interpretation of "T" as a doublon-hole continuum also explains that it is located at approximately U below the barycenter of the three-hole continuum. The latter is given by ω 0 = −6t for n = 2. With decreasing filling the support of the three-hole continuum shrinks, and so does its barycenter, but at a faster rate than the twohole continuum, so that the doublon-hole continuum "T" eventually merges with the doublon satellite for fillings approaching half-filling.
A stable, repulsively bound triplon forms in the presence of an additional repulsive nearest-neighbor Coulomb interaction V . This situation is shown in Fig. 7. Starting from the doublon-hole continuum at U = 6t and V = 0, and switching on V leads to the emergence of two satellites splitting off for sufficiently strong V . Both have the same shape as the non-interacting density of states ρ 0 (x) and are separated by about 2V from the continuum. The presence of two satellites rather than one is due to the fact that a triplon has an internal degree of freedom: On two sites, the three-hole states (1, 0) T ≡ c † 1↑ |vac. and (0, 1) T ≡ c † 2↑ |vac. form two configurations. When coupled by the Hamiltonian they result in a pair of bonding/antibonding eigenstates at energies E ± = U + 2V ± t. Note that with H 2 , we count the interaction between holes rather than electrons. The energy splitting of 2t very well explains the energy difference between the positions of the two satellites in the spectra in Fig. 7. Their widths are different, but become equal in the strong-coupling limit U, V t and are then given by 8t 2 /(U + 2V ). Fig. 8 shows the emergence of the two triplon satellites in the two-hole spectrum at n = 1.8 as obtained by DMRG. With V > 0, we find a much stronger growth of the entanglement entropy during the Chebyshev iteration. This is to be expected, as the model becomes non-integrable. 56 These spectra have therefore been calculated at a moderate energy resolution δE = 0.4t. For V = 3, for example, the MPS bond dimension at this point already exceeds 6200.  For V = 0, we find the dominating doublon satellite around ω ≈ −10t and the two peaks of the doublon-hole continuum at higher binding energies, separated by about 16t 2 /U ≈ 2.67t. One is located at ω ≈ −13.1t, while to other one is next to the intense van Hove singularity of the doublon structure at ω ≈ −10.4t, and barely distinguishable from it. It requires a high-resolution calculation to uncover this feature (δE = 0.04t, the purple curve in Fig. 8). For V > 0, the double-peak structure of the triplon appears, with a separation of 2t and shifting by about 2V with increasing V , as in the case of the threehole spectrum (Fig. 7). The doublon satellite, on the other hand, shows a considerable loss of spectral weight, but has a basically unchanged position when increasing V . In addition, we observe another resonance emerging out of the band-like part at ω ≈ −4t − V , which we can identify as two neighboring holes repulsively bound by V as a dimer.
The available resolution does not allow us to determine the fine structure of the triplon peaks. The linearprediction technique, 30,43,44 which is commonly used to enhance the resolution, has not provided us with additional spectral information in this case. Instead, we go to an even higher filling, n = 1.98, where an exact calculation is still feasible for a fairly large system. The result for L = 100 and δE = 0.01 is shown in the inset of Fig.  8. We observe that while the first of the triplon peaks appears right away, the second one only rises above the continuum when V exceeds a critical value of V c ≈ 1.3. Note that the shape of the former markedly differs from the free density of states, indicating that in contrast to the three-hole spectral function (Fig. 7), the triplon in A 2−hole,σ (ω) does not become stable in the whole Brillouin zone.

F. The quadruplon
The two-hole spectrum of the half-filled Hubbard model (see Fig. 3, lowest panel) displays a broad peak "Q" centered around ω ≈ −7.4t.
A two-hole excitation requires a double occupancy at a site i in the ground state 0, N . For strong U and N = L, a corresponding configuration is generated by a virtual second-order hopping process, which at the same time favors antiferromagnetic spin correlations. This is the famous Anderson super-exchange. 57 Starting from this initial-state configuration, the resulting configuration contributing to the N − 2-electron final state n, N − 2 must therefore contain two nearestneighboring empty sites, i.e., two neighboring doublons (see the sketch given in Fig. 9). One may ask whether these bind to form a heavier "quadruplon".
Before addressing this question quantitatively, let us attempt a simple argument favoring an interpretation of "Q" as a bound quadruplon opposed to two independently propagating doublons: Close to half-filling the low-energy states of the system exhibit strong antiferromagnetic correlations. Hence, there are two missing magnetic bonds associated with a single hole propagating in this antiferromagnetic background corresponding to an excitation energy 2J. Consequently, letting the two neighboring doublons, i.e., the two neighboring holes in the final state (Fig. 9b) dissociate and propagate independently costs 4J while forming a compound quadruplon requires 3J only. In the strong J-limit of the t-J model, Eq. (22), where the loss of kinetic energy upon quadruplon formation can be neglected, this argument 35 is strict. For weak J, related to the present case of the Hubbard model, it serves as a simple picture for the binding mechanism only.
A simple inspection of the presence of a binding energy can be done using the position of "Q": A two-doublon continuum would be expected at twice the energetic center of the doublon satellite, an additional binding energy  would mean a shift from that position. Looking at the closest dataset to half-filling (n = 62/60 in Fig. 3), we do not see such a shift within plotting accuracy, the "Q" peak at ω ≈ 9.8 coincides with the expected continuum center. Going to higher fillings, however, we find a shift that increases roughly linearly with n (from ∆ω ≈ 0.1 at n = 1.2 to ∆ω ≈ 0.4 at n = 1.6) towards lower binding energies, i.e. the doublons seem to bind attractively as the filling increases.
For n → 2 the phase space for quadruplon formation shrinks to zero. However, we can address the quadruplon formation in the final state and carry out a similar analysis as in Sec. III E by calculating the following four-hole spectral function for the completely filled band. The result is shown in Fig. 10.
For U = 0, the spectrum consists of a four-hole continuum at low binding energies; its support is given by 4W = 16t. For strong U , this continuum has extremely small spectral weight. Some weight has been transferred to an extended structure at intermediate excitation energies around ω = −8t − U which we identify as the one-doublon-two-hole continuum. Most of the weight, however, is taken by the correlation satellite located at ω ≈ −8t − 2U which corresponds to final states with two doublons. Notably, its shape does not resemble the self-convolution of a doublon satellite (cf. Fig. 1), as one would expect for a doublon-doublon continuum. This is in fact corroborated by the analysis of an effective model 14,58 , valid for strong U in the subspace with a fixed number of doublons, but no singly occupied sites: The model describes inherently stable bosons with the creation operator d † i = c † i↑ c † i↓ and a hard-core constraint (d † i ) 2 = 0 enforcing that any site is at most singly occupied with a single boson. Recall that throughout the paper, a doublon refers to a pair of holes rather than particles, and hence a doublon at site i corresponds to a hole n d i = 0, where n d i = d † i d i is the boson occupationnumber operator. The effective model includes a hopping term with amplitude J/2 = 2t 2 /U and a nearestneighbor density-density interaction of strength J, thus twice the hopping amplitude, but only half the bandwidth 2J. Note that this effective doublon-doublon in-teraction is attractive.
Our four-hole problem now becomes a two-hole problem in the effective model, Eq. (27), and can be solved in much the same way as the two-hole problem in the Cini-Sawatzky case, see Ref. 59 and Sec. III A. Separating the relative from the center-of-mass motion, the eigenvalues can be calculated for rather large systems. One does indeed observe a split-off resonance for certain values of the wave vector k, corresponding to a bound quadruplon. However, due to the small interaction strength J, the bound state is energetically positioned within the doublon-doublon continuum.
In order to demonstrate this, rather than showing just the eigenvalues, we calculate the k-resolved two-doublon spectral function within the effective model: Note that it is for simplicity defined such that the lowest binding energy is ω = 0. Fig. 11 shows the result of a calculation for N d = L = 1000 and J = 2t/3 corresponding to filling n = 2 and interaction U = 6t in the original model. The two-doublon spectrum corresponds to the satellite around ω = −8t − 2U − 2J ≈ −21.3t in Fig. 10. The support is given by ∆ω = 8t/3 ≈ 2.67t, though with very little weight on the high-binding-energy side. This is the expected value for two independent doublons: ∆ω = 2 × 4 × J/2. The k-resolution uncovers that the continuum actually has very small spectral weight and that the spectrum is in fact dominated by the intense quadruplon peak. Around the Brillouin-zone center the quadruplon is a resonance with a finite width, but still its weight dominates the whole energy range. For |k| 3π/4, it is fully split off from the continuum and has zero intrinsic width (the width is due to the finite δE only). We conclude that the quadruplon excitation represents a stable compound object formed by two bound doublons when its wave vector is close to the zone boundary, while it is a well-defined resonance with finite lifetime around the zone center.
The stability of the quadruplon can be explained by a phase-space argument: For a state from the two-doublon continuum, the hard-core constraint is not effective in the thermodynamic limit. Hence, the total energy of two doublons with wave vectors p and k − p is obtained as J cos(p) + J cos(k − p). Right at k = ±π, this adds up to zero, so that a bound quadruplon of arbitrary energy ω = 0 and wave vector k = ±π cannot decay. Note that the quadruplon branches split off from the continuum on the low-binding-energy side, as the quadruplon is bound by an attractive interaction of strength J, see Eq. (27). Indeed, the quadruplon binding energy at k = ±π is given by J = 2t 2 /3U ≈ 0.67t. 59 The same physics is found when considering the k-resolved four-hole spectral function in the full model: This is shown in Fig. 12 for U = 6t. Due to the now allowed singly occupied sites and a finite U , the one-doublon-twohole continuum around ω ≈ ω = −8t − U acquires some spectral weight. This corresponds to final states where one of the doublons has decayed. Most of the spectral weight, however, is taken by the quadruplon resonance or bound state at higher binding energies. For the quadruplon binding energy at k = ±π, we obtain about 0.58t due to the moderate strength of U = 6t.
At this point, we briefly note that the quadruplon is stabilized in the entire Brillouin zone by a sufficiently strong nearest-neighbor Coulomb interaction V > 0, very similar to the case of the triplon excitation discussed in Sec. III E. Fig. 13 shows the (local) four-hole spectrum for U = 6t and different V . Here, the quadruplon satellite, with a shape approaching the shape of the free density of states for strong V , is fully split off from the (very weak) two-doublon continuum and shows up on its high-binding-energy side since the interaction is repulsive. With increasing V , the quadruplon peak shifts with the additional binding energy 4V . Next, we continue to tackle the question of whether there is a stable bound quadruplon in the two-hole spectrum. First, we note that the structure "Q" in A 2−hole (ω) (see Fig. 3) continuously connects to the quadruplon peak analyzed with A 4−hole (ω). This is demonstrated in Appendix B where the (local) four-hole spectrum is presented in the entire filling range 1 ≤ n ≤ 2.
Hence, the question to be answered is the following: Do the two doublons in the final states contributing to "Q" for fillings n < 2 bind and form a compound object, i.e., is there an enhanced probability to find the doublons at neighboring sites? Above we have already given an argument relying on the energy shift of "Q" with respect to the expected position of the two-doublon continuum. However, a more direct analysis should be possible as the numerical approach is based on the many-body wave function which, in principle, contains all information on the system. One should therefore be able to extract specific information on the final states within a given frequency range. An implementation of this idea within the Chebyshev expansion method requires a new technique as described below.

G. Filter-operator technique
To extract specific information on the states contributing to the two-hole spectrum at a given frequency ω, we proceed in the following way: First, we define a state Ψ(ω) as the sum over the complete set of energy eigenstates weighted by the matrix element of the transition operator: Ψ(ω) is constructed to represent a typical final state contributing to the spectrum at frequency ω. Note that it is coincidentally identical to the "correction vector", which is employed within dynamical DMRG 36 to calculate the spectral function via Here, we employ it in a different way, choosing a Hermitian (but otherwise in principle arbitrary) "filter operator" F and calculate the following expectation value: This constitutes a "filtered" two-hole spectral function which we call A 2−hole [F ] (ω). The additional weight factors m F n give us some characterization of the eigenstates involved, depending on the choice of F . The normal spectrum is obviously recovered for F = 1. This quantity is numerically accessible as the diagonal part (ω = x = y) of the two-dimensional Chebyshev expansion defined by where the moments are given by in terms of the Chebyshev polynomials of the rescaled Hamiltonian, see Sec. II.
To convince ourselves that this method works, we first consider the well-understood limit of the completely filled band n = 2. We define an operator h i which projects out 14: Two-hole spectral function for n = 2 (thick black line) and two-hole spectra (normalized to area) using different filter operators (colored lines), calculated with exact states for U = 6t, L = 1000, δE = 0.01t (M = 1121).
states with an empty site i as and an operator s i , projecting out states with single occupancy at i: With this we can use F = 1 L i h i to filter out states characterized by doublons, and F = 1 L i s i to filter out states characterized by the absence of doublons, but the presence of single holes.
Results for U = 6t are shown in Fig. 14. We do indeed observe that the spectral function obtained with the former filter has support at the doublon satellite only, while the support of the latter coincides with spectral range of the band-like part and of the satellite. That the spectral function filtering out final states with singly occupied sites has finite weight at the doublon satellite reflects that doublon propagation in fact involves intermediate configurations with singly occupied sites.
Next we test the triplon case. Recall that the two-hole spectral function is computed at the central site of the chain, i.e., we take i = L/2 for the transition operator c i↑ c i↓ . To look for triplons, we can thus choose local filters and take F = h L/2 s L/2+δ at different distances δ. An enhanced weight at δ = 1 would indicate that the doublon and the additional hole bind and form a stable triplon.
The result of calculations for n = 1.8 is shown in Fig.  15. We do indeed see that the case δ = 1 does not look any more special than δ = 2, δ = 3 or δ = L/2. All choices roughly reproduce the shape of the two-hole spectral function and are barely distinguishable in the range of the peak at ω ≈ −13.1t, earlier identified as the triplon structure. This corroborates our previous result that doublon and hole do not bind. On the other hand, a stable triplon is formed for a nearest-neighbor interaction V = 1. This is indeed supported by the filter-operator technique: For V = 1 (see inset of Fig. 15) the case δ = 1 is in fact very different as compared to δ = 1: The weight of the filtered spectrum is strongly enhanced at frequen- cies around ω = −14.4t, where the peak of the bound triplon is located (cf. Fig. 8).
Finally, to investigate the quadruplon, we consider the filling n = 62/60 right above half-filling. The quadruplon peak is located at ω ≈ −9.8t below the doublon structure with barycenter at ω ≈ −6t, see Fig. 3. We have computed filtered spectra using F = h L/2 h L/2+δ .
Results are shown in Fig. 16. In this case, a filter with δ = 1 produces a strong signal at the quadruplon position even for V = 0, while the spectra for larger distances δ = 1 are again close to each other and roughly coincide with the original spectral function. This is a clear indication that there is a bound quadruplon, the contribution of which dominates the two-hole spectrum A 2−hole (ω), as compared to the doublon-doublon continuum.
H. Infinite doublon lifetime for k = π Let us return to the doublon. As has been discussed in Sec. III C, a stable excitation with infinite lifetime is found in the simple n = 2 limit and explained as a repulsively bound pair of holes. This interpretation breaks down for fillings close to half-filling, where the two-hole excitation propagates through a background of singly occupied sites with antiferromagnetic correlations. We are thus left with the question whether there is a stable doublon for fillings n < 2, and close to half-filling in particular.
The U -dependence of the doublon lifetime has been addressed in a number of previous theoretical as well as experimental studies, 10-14,58,60 starting from the concept of repulsive binding. Here, we would like to emphasize that the k-resolution of the local two-hole excitation adds an important new view to this problem. For the case of periodic boundary conditions, the k-resolved two-hole spectral function is defined as: Exploiting translational symmetry, it is straightforward to show that k-summation yields Hence, with A 2−hole (ω, k), one actually decomposes the local two-hole excitation into a coherent linear superposition of k-dependent two-hole excitations: that the state d k 0, N has a strongly k-dependent lifetime. At the zone boundary k = π, the lifetime is even infinite. This can be verified numerically and understood analytically. Let us first discuss the numerical results. Eq. (37) can be used to approximate the two-hole spectrum with a calculation for a system with open boundaries. In this case, the k-summation of A 2−hole (ω, k) [Eq. (38)] yields the local spectrum, but averaged over all sites, rather than the local spectral function at the central site i, as defined in Eq. (5). For a system with size L = 60, however, the approximation is already excellent. Note that this can also be thought of as carrying out a Fourier transform with two momenta, A 2−hole (ω, k, k ), but then only looking at the diagonal component with k = k : A 2−hole (ω, k) ≡ A 2−hole (ω, k, k). Fig. 17 shows the k-resolved spectrum for n = 1.8 and U = 6t. The band-like part extends from ω = 0 to ω ≈ −8t with most of the spectral weight at low binding energies and a ridge-like structure around k = 0. Its overall spectral weight, however, is very small compared to the weight of the doublon satellite, which is centered around ω ≈ −10t. Interestingly, one finds that the doublon satellite is not broadened at all for k = ±π. In fact, the structure at k = ±π should be interpreted as a δpeak at frequency ω = ω 0 (π) ≈ −9.25t. Note that on the logarithmic scale, there are seemingly sidebands visible close to ω 0 (π) and k = ±π. These are, however, simply the numerical artifact of Gibbs oscillations. As shown below, the position of the δ-peak is precisely given by ω 0 (π) = U − 2µ ≈ −9.25t, where the value of µ is obtained from the ground-state calculation [Eq. (7)]. When |k| < π, the satellite acquires some finite width, i.e., the doublon transmutes from a stable bound object into a resonance. Its lifetime is shortest around k = 0. This is opposed to the simple n = 2 limit, where the doublon is stable in the entire Brillouin zone. The situation at a filling n = 1.1 is very similar. Results are shown Fig. 18. Here, the k-resolved two-hole spectrum consists of an intense doublon satellite (centered around ω ≈ −5t) and the quadruplon structure (around ω ≈ −10t), as well as a band-like part with nonzero, but extremely low weight close to ω = 0 which is not visible in the figure. Referring to the results obtained by means of the filter-operator technique in Sec. III G, we interpret the quadruplon structure as an object composed of two doublons but with a finite lifetime, i.e., a resonance. Most of its weight is located close to k = 0. It is tempting to assume that the quadruplon is entirely stable at k = π, similar to the four-hole spectrum case (Fig.  12), but the whole spectral weight at the zone boundary is instead taken by the doublon satellite. We have already seen a similar effect in the case of the triplon: In the three-hole spectral function for n = 2 (Fig. 7), it clearly has a different dispersion from the two-hole spectral function for n < 2 (the inset of Fig. 8), while the energetic positions follow the same rules.
Turning to the doublon in Fig. 18, we observe that it has infinite lifetime at k = ±π again, but now also displays a distinct "shadow" around k = 0. For n = 1.1 we find 2µ = 9.25, and thus expect the doublon at energy ω 0 (π) = U − 2µ ≈ −3.25t. This perfectly fits with the position that can be read off from Fig. 18.
The stability of the doublon at the zone edges can be understood analytically by referring to the "hidden" charge-SU(2) symmetry of the Hubbard model. 61,62 We define For the model considered here, and in general for a Hubbard model with V = 0 on a bipartite lattice in D dimensions, one easily verifies that η is an eigenoperator of H and of the total particle-number operatorN : η is related to the total isospin T , the components of which are given by T x = (η + η † )/2, T y = (η − η † )/(2i) and T z = (L −N )/2, and represent the generators of a global charge-SU(2) symmetry of the model at halffilling: [T , H − µN ] = 0 for µ = U/2. Starting from the total spin S, which generates the standard spin-SU(2) group, the isospin is obtained by a spin-asymmetric and staggered particle-hole (Shiba) transformation U sh with U † sh c i↑ U sh = (−1) i c † i↑ . We have U † sh SU sh = T . At half-filling (N = L) the ground state 0, N is a nondegenerate total spin singlet. This implies immediately that 0, N is a total isospin singlet as well. Hence, we have η 0, N = L = 0. Since d k=π = η/L for k = π, the transition operator in the Lehmann representation of the k-resolved two-hole spectral function, Eq. (37), is just proportional to η. For the case of half-filling this implies for all excitation energies ω. This is consistent with the DMRG calculations which predict the doublon satellite to vanish for n = 1, even in the entire Brillouin zone (see Fig. 3).
Off half-filling, the chemical potential µ = U/2 explicitly breaks the charge-SU(2) symmetry and produces a massive collective mode showing up as a single δ-peak in the two-hole spectrum at wave vector k = π (the same holds for k = −π). This can be inferred from Eq. (41). Namely, it is easy to see that, for arbitrary N , the state η 0, N is an exact eigenstate of H with energy E (N ) 0 −U . Inserting this into Eq. (37), yields This explains the observed position ω 0 (π) = U − 2µ.
Integrating the k-resolved spectrum, Eq. (37), over ω, yields the total weight at a given k: This gives k α(k) = i n i↑ n i↓ , consistent with the global sum rule Eq. (18). At k = π, as expressed by Eq. (43), the whole spectral weight condenses to α(π) at frequency ω = U − 2µ. The weight of the mode at the zone boundary is then given by α = L −1 0, N η † η 0, N , see Eqs. (40) and (44), and can be expressed in terms of the isospin as follows: where we also used M T = (L − N )/2. We conclude that the weight of the mode linearly increases from zero at half-filling to unity at n = 2. This explains the strong difference in the weights seen in the DMRG results (Figs. 17 and 18). The same analysis can be done for the k-resolved twoparticle spectrum and fillings n > 1. Here, however, the exact eigenstate η † 0, N would lead to a δ-peak at negative frequency ω = ω 0 (π) = U − 2µ < 0 for n > 1 [which would contradict Eqs. (6) and (7)] but the weight α = L −1 ηη † = L −1 T 2 − T z (T z − 1) vanishes exactly since M T = −T . For fillings n < 1, the δ-peak in the two-particle spectrum shows up on the ω > 0 side with a non-zero weight since M T = +T in this caseanalogous to the two-hole spectrum for n > 1 and related to the latter by particle-hole transformation.
Physically, the collective mode describes the coherent propagation of a doublon with energy U − 2µ through the lattice without any scattering, i.e., the final state with wave vector k = π in the two-hole spectroscopy has infinite lifetime. This is enforced by the remaining isospin-rotation symmetry around the z-axis which protects the "η-mode" at the zone edge from decay. Note that the η-mode only shows up off half-filling and is not bound to off-diagonal long-range order as has been discussed previously. 61,62 For bipartite lattices of higher dimensions, the η-mode is found at the Brillouin-zone edges k = (±π, ±π, ...).

IV. DISCUSSION
While the spectrum of one-particle excitations of the Hubbard model has been discussed extensively over the recent decades, much less is known for excitations involving two or more particles or holes. Referring to the Lehmann representation of the spectral densities, it is obvious that the different spectroscopies are intimately related since they involve the same set of final states (when starting from a ground state with properly chosen total particle number). The main difference, however, is in the spectral weight attached to the final states. This strongly depends on the type of excitations studied. Dif-ferent spectroscopies therefore provide largely different perspectives on the same physics.
For strong U , the one-hole spectrum of the Hubbard model (corresponding to photoemission) exhibits the famous Hubbard bands. In the high-density regime (filling close to n = 2), the lower Hubbard band is built from final states involving configurations with two holes at the same site. It is therefore shifted by about U to higher binding energies and provides an interpretation of the Ni-6eV satellite, 63 for example. The weight of the satellite is small since the probability to excite the photoelectron from a singly occupied site is low for n → 2.
The same final states have a much higher weight in the two-hole spectrum (corresponding to Auger-electron spectroscopy) and the related peak is referred to as the Cini-Sawatzky satellite. This satellite in the two-hole spectrum can even be studied in the limit of the completely filled band (n = 2) where, in the one-hole spectroscopy, the LHB disappears or has zero weight. The related two-hole problem has attracted a lot of interest as it can be solved analytically and provides the concept of a repulsively bound pair of holes (called "doublon" in this paper).
Besides recent experiments addressing doublon physics in systems of ultracold atoms trapped in optical lattices, the two-hole problem has a longer history and has in particular stimulated much work on high-resolution Auger spectroscopy from metals. For n < 2, however, a comprehensive overview of the two-hole excitations, or even an understanding of the spectral features has not yet been achieved. The obvious reason is that for n < 2, one is confronted with an intractable many-body problem in general.
For moderately large one-dimensional lattices, however, one can profit from the recent advances of numerical techniques based on matrix-product states. As is demonstrated here, this allows to study the two-hole excitation spectrum to a high degree of accuracy. We have employed a matrix-product state implementation of the density-matrix renormalization group in combination with the Chebyshev polynomial expansion technique to compute different excitation spectra. As the excitation is local (k-resolved spectra are computed by combining different local excitation spectra) the growth of the entanglement entropy with expansion order scales logarithmically and is very moderate (though it is larger for higher-order spectroscopies and quite substantial for a nearest-neighbor interaction V = 0). Going to high orders (typically M ∼ 10 3 ) provides us with the numerically exact solution on a frequency scale larger than a "resolution" δE (typically δE ∼ 10 −1 t).
Computing the two-hole spectrum of the Hubbard model with this technique uncovers a couple of new phenomena. Some are intrinsically related to the onedimensionality of the lattice, such as spin-charge separation. These will be discussed in an upcoming paper but have been disregarded here as far as possible. With this we concentrate on those aspects which should, at least qualitatively, be the same for higher dimensions as wellwith the notable exception of physics related to ground states with spontaneously broken symmetries.
We find that the overall two-hole spectrum is governed by the physics of multiplons, i.e., bound states of infinite lifetime or resonances of finite spectral width showing up at characteristic excitation energies and with characteristic dispersions and weights. The footprints of multiplons are actually ubiquitous in the one-, two-and moreparticle or hole spectra of the (extended) Hubbard model. This multiplon physics is clearly developed at strong U (U = 6t = 1.5W has been chosen throughout the paper) where the typical energy scales are well separated. In this regime some charge fluctuations can be disregarded, and effective models, like an effective hard-core-boson model (at high fillings) or the t-J model (fillings close to halffilling), are very useful.
For n = 2 and strong U , most of the weight of the two-hole spectrum is concentrated in the doublon peak. The two additional holes in the final state form a compound object of infinite lifetime at binding energies approximately centered around U below the barycenter of the two-hole continuum. For fillings n → 1, however, the interpretation of the doublon peak as a repulsively bound hole pair changes, and it rather reflects the motion of an empty site through a lattice of singly occupied sites with antiferromagnetically correlated magnetic moments, as described by the t-J model.
While the doublon is well defined in the whole filling range at U = 6t, it cannot be expected to be stable (for n < 2) and should decay by emitting particle-hole excitations. With decreasing n, the phase space for this decay channel increases; this destabilizes the doublon. With increasing U , energy conservation requires the emission of multiple particle-hole excitations in a high-order scattering event; this stabilizes the doublon. Besides these overall trends, however, the doublon lifetime has turned out to be strongly k-dependent. Interestingly, it is even infinite at k = ±π -for any filling n > 1 and for bipartite lattices of arbitrary dimension. This phenomenon is related to the "hidden" charge-SU(2) symmetry of the Hubbard model at half-filling. Off half-filling, the coupling of the chemical potential µ to the z-component of the isospin explicitly breaks this symmetry resulting in a massive collective mode in the two-hole spectrum (the doublon at k = ±π) with energy ω = U − 2µ and with weight α = n − 1.
The two-hole spectrum also uncovers heavier multiplons which show up in certain filling ranges with lower weight. Besides the doublon, we have identified a triplon and a quadruplon excitation. The triplon consists of a hole pair (a doublon) with an additional hole on a nearest-neighboring site. The quadruplon is a pair of neighboring doublons. Let us note in passing that those can even be found, with extremely low weight at U = 6t, in the one-hole spectrum. Furthermore, the multiplon hierarchy can also be continued: There is for example evidence for a quintuplon, i.e., five holes on three sites, showing up as a weak peak on the high-binding-energy side of the quadruplon in the four-hole spectrum (see Appendix B).
For a Hubbard-type model with local interaction and with two (spin) degrees of freedom per site, an odd-even effect is a natural expectation. It would seem that multiplons composed of an odd number of holes do neither form stable compounds nor well-defined resonances unless a nearest-neighbor density interaction V is added to the Hamiltonian. Contrary, the doublon and the quadruplon appear as resonances with comparably small width or even as stable objects with infinite lifetime in the kresolved two-hole spectrum. The four-hole spectrum at n = 2 represents an interesting special case, where a stable quadruplon is found in a finite k-range around k = ±π until it is immersed in the doublon-doublon continuum and becomes a sharp resonance of finite width. Odd multiplons, when stabilized by V , possess one or more internal degrees of freedom as there are hopping processes which do not affect their integrity in case of strong interactions. In the case of the triplon, for example, this manifests itself in a double-peak substructure on a scale 2t.
It is important to distinguish clearly between dispersion and decay. The width of multiplon peaks in the local two-hole spectrum is mainly due to dispersion, which is finite for finite U < ∞ and can be understood best in the framework of effective models obtained from strongcoupling perturbation theory. Also, k-resolution helps to identify the intrinsic lifetime broadening of resonances.
To distinguish between resonances and continua can be problematic in some cases. Here, the triplon and the quadruplon represent good examples, as their excitation energies lie within the doublon-hole or within the doublon-doublon continuum, respectively. However, for V = 0 only the quadruplon is a compound object with finite lifetime, while the triplon decays immediately. To analyze such cases and come to unambiguous conclusions, we have applied a new filter-operator technique. This allows to extract information on the states contributing to the spectrum at a given frequency, which is specifically dependent on the choice of an arbitrary filter operator. Within the framework of Chebyshev polynomial expansion, this can be implemented in a straightforward way.
There are several lines along which one can extend the present study, as for example study the manifestation of spin-charge separation in the two-hole spectrum, or the footprints of spontaneous U(1)-symmetry breaking for U < 0. Another set of questions regards the real-time dynamics of multiplons and the relation to time-dependent (pump-probe) Auger-electron spectroscopy. Work in these directions is in progress.

acknowledged.
Appendix A: (Inverse) photoemission spectra Fig. 19 shows the (local) one-hole excitation spectrum [see Eq. (23)] and the one-particle excitation spectrum, at U = 6t and for fillings ranging from half-filling (n = 1) up to limit of the completely filled band (n = 2). These correspond to photoemission (one-hole spectrum) and inverse photoemission (one-particle spectrum). For n = 2, the unrenormalized density of states is recovered: A 1−particle,σ (ω) = ρ 0 (ω + µ). The van Hove singularities at ω = 0 and ω = −4t are somewhat broadened due to the finite δE = 0.2t. For fillings n < 2, the spectrum consists of the upper Hubbard band (UHB) at low binding energies and the lower Hubbard band (LHB) showing up at higher binding energies, approximately at U below the UHB. Note the strong filling-dependent spectral-weight transfer. At half-filling the LHB and the UHB are placed symmetrically to ω = 0. The substructure of the LHB and UHB can be analyzed and understood in terms of collective spinon and holon excitations, see Refs. 38-40. Appendix B: Four-hole spectra Fig. 20 shows the (local) four-hole excitation spectrum at U = 6t for various fillings in the range of excitation energies ω where the structure "Q" is seen in the two-hole spectrum (Fig. 3). The position of the main peak in Fig.  20 approximately coincides with the position of "Q" in the two-hole spectrum for all fillings n, if the same final subspaces are compared (for example, at ω ≈ −14.5t in A 2−hole (ω) for n = 84/60 = 1.4 compared to ω ≈ −14.6t in A 4−hole (ω) for n = 86/60 ≈ 1.43). For n → 2, the main peak evolves into the quadruplon excitation that has been analyzed in Sec. III F. The weight of the quadruplon strongly depends on the spectroscopy considered. In the two-hole spectrum, it is most intense for n → 1, while it has maximum weight for n → 2 in case of the four-hole spectroscopy.
The latter also uncovers the existence of a "quintuplon", i.e., a final state composed of two neighboring doublons plus an additional neighboring hole. This is visible as the additional structure on the high-bindingenergy side of the quadruplon in Fig. 20 (see ω ≈ −22t for n = 1.8, for example). As in the triplon case, we suspect that the quintuplon can be stabilized as a truly bound state with infinite lifetime when switching on a nearestneighbor Coulomb interaction. For V = 0, however, it is not stable and should correspond to a hole-quadruplon continuum of final states.