Finite-temperature photoemission in the extended Falicov-Kimball model: a case study for Ta$_2$NiSe$_5$

Utilizing the unbiased time-dependent density-matrix renormalization group technique, we examine the photoemission spectra in the extended Falicov-Kimball model at zero and finite temperatures, particularly with regard to the excitonic insulator state most likely observed in the quasi-one-dimensional material Ta$_2$NiSe$_5$. Working with infinite boundary conditions, we are able to simulate all dynamical correlation functions directly in the thermodynamic limit. For model parameters best suited for Ta$_2$NiSe$_5$ the photoemission spectra show a weak but clearly visible two-peak structure, around the Fermi momenta $k\simeq\pm k_{\rm F}$, which suggests that Ta$_2$NiSe$_5$ develops an excitonic insulator of BCS-like type. At higher temperatures, the leakage of the conduction-electron band beyond the Fermi energy becomes distinct, which provides a possible explanation for the bare non-interacting band structure seen in time- and angle-resolved photoemission spectroscopy experiments.


Introduction
In solid state physics, superconductivity is a classic example of emerging quantum coherence on a macroscopic scale. An essential prerequisite for this is the pairing of electrons, normally of opposite spin and close to Fermi surface of metals. Similarly, valence-band holes and conducting-band electrons may form excitonic pairs in semimetals with small band overlap or in semiconductors with a small band gap, triggered by their Coulomb attraction, where the excitons are expected to condense in a macroscopic coherent state at low temperatures under restrictive conditions [1][2][3][4][5]. This state is called "excitonic insulator" (EI) because it exhibits no "super-transport" properties. It is stabilized by the opening of a gap in semimetals or by the flattening of the valance-band top (conduction-band minimum) in semiconductors, which both can be detected by angle-resolved photoemission spectroscopy (ARPES). Usually exciton condensation in a semiconductor setting is discussed in terms of Bose-Einstein condensation (BEC), while those appearing for semimetallic (noninteracting) band structure should be described in close analogy with the Bardeen-Cooper-Schrieffer (BCS) theory [4][5][6][7][8][9].
Even though the obvious EI scenario was predicted almost 60 years ago, it was only recently that some material classes have shown some signatures of an EI ground state. In this respect one of the most promising candidates seems to be the quasi-one-dimensional material Ta 2 NiSe 5 , which possesses a tiny direct gap at the Γ point of the Brillouin zone. At a critical temperature T c ≈ 328 K, Ta 2 NiSe 5 encounters a second-order transition from an orthorhombic (T > T c ) to a monoclinic (T < T c ) structure [10]. ARPES experiments revealed a flattening of the valence band when cooling below the transition temperature [11,12], which has been taken as indication of an EI state within a model simulation [13]. Motivated by this, various experiments have recently performed for this material, using NMR [14], inelastic x-ray scattering [15], scanning tunneling microscopy or spectroscopy [16,17], and time-and angleresolved photoemission spectroscopy (t-ARPES) [18][19][20]. Besides there is the isostructual compound Ta 2 NiS 5 , which opens up the possibility to perform comparative measurements and analysis, not only of the photoemission [21] but also of the magnetic susceptibility [10] and the optical conductivity [22,23].
Unfortunately, some experimental findings and theoretical predictions are inconsistent so far. For instance, the optical conductivity spectra in Ta 2 NiSe 5 [23] show an extra peak in the low-energy regime (at about ω 0.4 eV) that cannot be reproduced by the densityfunctional-theory-based simulation [24]. The authors of Ref. [23] attributed this low-energy peak to the giant oscillator strength of spatially extended exciton-phonon bound states, while the theoretical work [24] concluded that the interorbital Coulomb interaction between valence and conduction bands should be sufficient to explain this peak. Another important issue is whether the formation of the EI in Ta 2 NiSe 5 follows a BCS or a BEC scenario. While the small value of the transport gap [22] implies that the excitonic state in Ta 2 NiSe 5 is of BCS type, the spectral weight broadening of the flat band with increasing temperature has been taken as a signature of a breaking of the quasiparticle peak structure, i.e., of dissolving a Bose-Einstein condensate [13].
In contrast to the advanced experimental investigations performed for Ta 2 NiSe 5 in the last few years, most of the theoretical studies are mainly based on weak-coupling approaches, especially if the single-particle spectra have been addressed [13,25,26]. When focusing on strictly one-dimensional (1D) systems this is not necessary, however, because we can calculate both ground-state and dynamical quantities for the considered model Hamiltonians by rather unbiased numerical techniques like the density-matrix renormalization group method (DMRG) [27] and its further developments. For sure, such a de-facto approximation-free treatment is essential to resolve some of the issues carved out above.
To this end, in this work, we re-examine the properties of the ground state of the 1D extended Falicov-Kimball model (EFKM), as well as the behavior of various spectral functions at zero and finite temperatures, by means of the time-dependent DMRG (t-DMRG) technique [28,29] in the matrix-product-state (MPS) representation. The EFKM at half band filling can be considered as the perhaps minimal model for Ta 2 NiSe 5 . We demonstrate that the best-suited parameter set typifies Ta 2 NiSe 5 as an EI in the BCS regime; here, the zerotemperature photoemission spectra exhibit two-peak structure. At finite temperatures the gap melting is clearly visible in the photoemission spectra, and the leakage of the conduction band beyond the Fermi energy becomes distinct as temperature increases, which might be related to the bare band structure obtained in the t-ARPES experiment [19].
The paper is organized as follows. In section 2, we introduce the EFKM Hamiltonian, present its ground-state phase diagram, including the BCS-EI BEC-EI crossover regime deduced from pair-condensation amplitude, and comment on the best model parameters used for describing Ta 2 NiSe 5 . Section 3 contains our numerical data for the photoemission spectra in the EFKM, which are discussed with regard to the experimental results for Ta 2 NiSe 5 . Our main conclusions are presented in section 3. Details of the used numerical scheme can be found in Appendix A. For comparison, Appendix B provides results for the photoemission spectra in the 1D half-filled Hubbard model.

Theoretical modeling 2.1 Extended Falicov-Kimball model
Let us first introduce the 1D EFKM [30][31][32][33], which has been argued to be the minimal theoretical model for Ta 2 NiSe 5 [26,34]. The Hamiltonian readŝ whereα † j (α j ) denotes the creation (annihilation) operator of a spinless fermion in the α = {c, f } orbital at Wannier site j,n α j =α † jαj , U is the local Coulomb repulsion between c and f electrons staying at the same lattice site, D parametrizes the level splitting of c and f orbitals, and µ is the chemical potential. Representing the orbital flavour of the EFKM by a pseudospin variable [30,35], t c → t ↑ t f → t ↓ ,ĉ j →ĉ j,↑ , andf j →ĉ j,↓ , the model (1) can be viewed as asymmetric Hubbard model with spin-dependent hopping in a magnetic field.  Figure 1: (a) DMRG ground-state phase diagram of the 1D half-filled EFKM. Star ( ) and diamond ( ) symbols mark the parameter sets that will be used in simulations of the spectral functions A(k, ω) in Sec. 3. (b) and (c): Condensation amplitudes F (k) at various n c j for U/t c = 1 and 8. Data obtained by DMRG with L = 60 and periodic boundary conditions. Dashed lines give the corresponding Fermi momenta k F = π n c j in the noninteracting limit. In all cases, t f /t c = −0.5.
In other words, the standard Hubbard model (see. Eq. (8) in Appendix B.1) is obtained by setting t c = t f and D = 0 in the EFKM. In what follows, we take t c as the unit of energy and, with a view to Ta 2 NiSe 5 , consider the half-filled band case only. Note that direct f -c hopping is prohibited by symmetry for this material [26].
The DMRG ground-state phase diagram of the half-filled 1D EFKM has been worked out previously [34]. Depending on the strength of the orbital level splitting, the system realizes one of three insulating phases: a state characterized by staggered orbital order (SOO) with n c j +n c j+1 = n f j +n f j+1 = 1, the EI where the c-f electron coherence emerges (here, 0 < n c j < 1/2 < n f j < 1), and the (rather trivial) band insulator (BI) where n c j = 0 and n f j = 1. While the EI-BI phase boundary is given analytically [36], the SOO-EI quantum phase transition has to be determined numerically, e.g., from where E 0 (N f , N c ) is the ground-state energy for a finite system with L lattice sites and N f f -electrons and N c c-electrons. Figure 1(a) shows the ground-state phase diagram of the 1D EFKM in the U -D plane for t f /t c = −0.5. Obviously, the SOO phase takes up only a small space in the vicinity of D = 0. Its boundary (3) can be easily extracted from a DMRG calculation of the various ground-state energies at fixed system size. The EI and BI phases are separated by the transition line (2).
As demonstrated in Ref. [34], the EI is characterized, by the behavior of the binding energy of electron-hole pairs [37] , and the condensation amplitude in the momentum space. In (4), |ψ 0 is the ground state |N f , N c with fixed numbers of fand c-orbital electrons, and |ψ 1 is the excited state |N f − 1, N c + 1 . If the electron-hole pairs are only weakly bound (E B /t c 1, BCS-EI regime), F (k) develops a sharp peak at the Fermi momentum k F = π n c j [cf. Fig. 1(b) for U/t c = 1]. Here, Fermi "surface" effects play an important role. On the other hand, for tightly bound excitons (E B /t c 1, BEC-EI regime), F (k) has a (broad) maximum at k = 0, see F (k) for n c j = 0.1 and 0.2 at U/t c = 8 in Fig. 1(c). One can therefore define a BCS-BEC crossover region, where F (k) has a maximum at 0 < k < k F . Within the the EI, the respective defined BCS-BEC crossover range is visualized in Fig. 1(a) by the shaded region. We see that for the t f /t c = −0.5 ratio used, an EI of BEC type appears for U/t c 3 in the vicinity of the EI-BI transition line only. In this respect we note that for U/t c = 1 [U/t c = 8] the EI-BI transition occurs at D c 2 /t c 2.16 [D c 2 /t c 0.54], and the results depicted in Fig. 1

Ta 2 NiSe 5 model parameters
To determine the optimal parameter set for an (extended) Falicov-Kimball-model-based description of Ta 2 NiSe 5 , ARPES data is mainly used. For this, Seki et al. [13] utilized a temperature-dependent variational cluster approach for the three-dimensional EFKM, assuming t c = −t f , however, not least because of the reduced computational costs due to particle-hole symmetry in this case. A perhaps more realistic approach was put forward by Kaneko et al. [26] who considered a three-chain electron-phonon-coupled system and exploited band-structure calculations together with a mean-field analysis. Here, the estimated parameter set is t c = 0.8, t f = −0.4, U 0.55 and D = 0.2 (in units of eV). In all these efforts, it is worth observing that band flattening detected in the ARPES experiments only occurs in a relative narrow region of momentum space, specifically for |k x | 0.1Å −1 (taking into account that the lattice constant of the chain direction is a = 3.51232Å [38], this means for |k|/a 0.35). Since for t f /t c = −0.5 and U/t c 1 the system is in the EI-BCS regime where the maxima of F (k) are almost equal to the peak position of the single-particle spectral functions [34], the flat band can appear only for |k| k F . This also implies that the strength of the level splitting should be very close to the EI-BI phase transition line, as for n c j = 0.1 where k F = 0.1π(< 0.35). Therefore, in this study, we discuss the spectral functions of the 1D EFKM using the parameter set t f /t c = −0.5, U/t c = 1 (which is slightly larger than above estimation U/t c = 0.55/0.8 0.69), and D/t c = 2.11. This allows us to carry out the DMRG simulations for a mean c-electron density n c j = 1/10 at T = 0. This parameter set is marked by a star ( ) in Fig. 1(a).
To provide a contrasting perspective, we will also consider the strong-coupling regime, specifically for U/t c = 8, where the BEC-type EI is realized in a wider region of the density (i.e., for n c j 0.25). Here, we use D/t c = 0.53, and again n c j = 1/10, for the calculation of the spectral functions at T = 0 [cf. the symbol in Fig. 1(a)].
Let us finally emphasize that we will approximate the doubly-degenerate conduction celectron bands [26] by a single band and neglect any electron-phonon coupling effects [39] on the EI formation.
3 Numerical results for the extended Falicov-Kimball model 3

.1 Spectral functions
The full single-particle spectrum A(k, ω) consists of two parts, where A − (k, ω) and A + (k, ω) denote the photoemission spectrum (PES) and inverse PES (IPES), respectively. In the case of the EFKM, A ± (k, ω) splits into contributions from the c and f electrons, i.e., A ± (k, ω) = α=c,f A ± α (k, ω). As a result, to determine A(k, ω) for t c = |t f |, we have to compute four dynamical correlation functions, A ± α (x, t) with α ∈ {c, f }. These circumstances distinguish the EFKM from the simpler Hubbard-model case, where only one of the four correlation functions A ± σ (x, t) with σ ∈ {↑, ↓} is needed to obtain full spectra because of particle-hole and spin-flip symmetries. To obtain A ± (k, ω) numerically, we simulate the dynamic correlation function in real space and then carry out a Fourier transform where · β denotes the expectation value at inverse temperature β. In doing so we exploit the t-DMRG technique together with the purification method [40,41] for MPS and so-called infinite boundary conditions (IBC) [42]. Moreover, within t-DMRG framework, we use a second-order Trotter-Suzuki decomposition with a time step τ = τ β = 0.01 for the real-and imaginary-time evolutions and a maximum bond dimension of 800, so that the truncation error will be smaller than 10 −6 in all the calculations presented. Appendix A gives a more detailed description of our numerical scheme, whose accuracy is tested for the Hubbard model in Appendix B.

Single-particle spectra at zero temperature
We now discuss the single-particle spectrum of the EFKM at half-filling for both zero (T = 0) and finite (T = 1/β) temperatures with a particular focus on the flattening of the band structure observed in ARPES experiments for the quasi-1D potential EI material Ta 2 NiSe 5 around the Γ point. Clearly, in a strictly 1D model system, a symmetry-broken EI phase can only exist at zero temperature. Nevertheless, the spectral properties of the 1D EFKM at low temperatures should reflect the existence and nature of an EI phase at T = 0. Figure 2 displays the zero-temperature, wave-vector and energy resolved spectral functions in the EFKM, using the parameter set best-suited for Ta 2 NiSe 5 : t f /t c = −0.5, U/t c = 1 and D/t c = 2.11 [marked by the star in Fig. 1(a)]. Here, panels (a) and (b) give the c-and f -orbital parts of the spectral functions A c (k, , respectively, which are combined in panel (c). In this regime, starting , and U/t c = 1. The mean celectron density is fixed to be n c j = 0.1. The dashed lines in panel (d) mark the momentum interval considered in ARPES experiments on Ta 2 NiSe 5 [11]. Results are obtained by t-DMRG using IBC and a window size L W = 200 to keep the values of the dynamical correlation functions on the boundaries less than 10 −8 up to the target time t fin · t c = 16. We note that the unit cell for the iDMRG simulations needs to be large (N uc = 10) to realize the small c-electron density n c j = 0.1.
out from a semi-metallic bare band situation, the EFKM realizes an EI phase of BCS type built up from weakly-bound electron-hole pairs. Since the Coulomb attraction is relatively weak, A c (k, ω) and A f (k, ω) almost follow the unrenormalized c-and f -band dispersions. Concomitantly, we find a more or less uniform distribution of the spectral weight and weak incoherent contributions close to band edge only. Nevertheless, a gap develops at the "bare" Fermi momentum k F = π n c j , which is exponentially small, however, and therefore difficult to detect in panels (a)-(c) in view of the used Lorentzian broadening η L /t c = 0.1. It is actually only confirmed in the zoomed-in panel for the photoemission part [ Fig. 2(d)]. Most notably the PES shows an "M"-shaped band structure with two peaks close to ±k F , just as detected in the Ta 2 NiSe 5 -ARPES experiments at low temperatures in the vicinity of the Γ point [13,20], even though the energy minimum at k = 0 turns out to be somewhat too low in our EFKM simulations. In this respect, better agreement with experiments might be reached by reducing the density of c-electrons further ( n c j < 0.1). In any case, the "flat-band" region falls entirely into the restricted momentum window monitored within the ARPES measurements, see dashed lines in Fig. 2(d).
If the attraction between electrons and holes is strong, tightly bound excitons will be formed. At the same time the Hartree shift enlarges the orbital splitting and drives the system further into a semiconducting situation [33]. As a result the EI (low-temperature) phase appears as BEC of preformed (excitonic) pairs [6]. This needs to be reflected in the spectral properties of the EFKM as well, which can not longer be described adequately by weak-coupling approaches, such as Hartree-Fock [43], random phase approximation [35] or (projector-based) renormalization methods [25]. Figure 3 shows the different single-particle spectra A c,f (k, ω) and A (−) (k, ω) for the strongcoupling model parameters belonging to the diamond symbol ( ) in Fig. 1(a). Quite clearly, now the essential characteristics of the spectra are a huge gap for single-particle excitations In view of the pronounced band flattening, the obvious question is whether Ta 2 NiSe 5 fits into such a EI BEC-type scenario, especially because this material is in truth semiconducting. The answer, however, is no, as can be seen from the results for A(k, ω) alone. In the BEC-EI regime of the EFKM, a two-branch structure emerges in the PES A − (k, ω) in the vicinity of the Brillouin zone center (Γ point), just as in the strong-coupling regime of Hubbard model (see Fig. 8(d) in Appendix B.2), which is not observed in ARPES. Furthermore, the flat-band region in momentum space appearing in the EFKM PES spectrum is notably wider than that for Ta 2 NiSe 5 ; cf., the zoomed-in panel Fig. 3(d).
To draw an interim balance, it seems that the ARPES experiments on Ta 2 NiSe 5 can be explained in the framework of the pure 1D EFKM with parameters deep in the weak-coupling region and very close to the EI-BI transition, suggesting an EI state of BCS type.

Photoemission spectra at finite temperatures
To further substantiate that the possible EI state of Ta 2 NiSe 5 might be of BCS type, we now analyze the single-particle spectral functions at finite temperatures. It is of specific interest to prove the leakage of spectral weight of the lower band to energies above the Fermi energy in the EFKM, similar to the melting of the Mott gap in the Hubbard model, see Ref. [44] and Appendix B.2. Such leakage reflects the intensity tail above E F after pumping observed in a recent t-ARPES experiment for Ta 2 NiSe 5 [19]. Likewise we should verify in terms of the EFKM that the upper band in the PES at T > 0 is missing or merged into a single band due to the small value of U , just as indicated by ARPES on Ta 2 NiSe 5 [11,12,21].
To achieve this technically, we perform finite-temperature simulations, working with a grand canonical ensemble and employing t-DMRG with IBC and an N uc = 2 (2 physical and 2 auxiliary sites) unit cell for the iMPS. The chemical potential µ is fixed by the infinite TEBD so as to fulfil the condition n c j + n f j = 1 for each target temperature, before starting  simulations for the dynamical correlation functions. Figure 4 provides the PES A − (k, ω) spectrum of the 1D EFKM for the same parameters as in Fig. 2] but at finite temperatures. Clearly, at low temperatures (T /t c = 0.1 [ Fig. 4(a)]) the results differ marginally from their zero-temperature counterparts given in Fig. 2]; only the peak height is slightly reduced. Because of the small charge gap, the gap melting occurs already at very low temperatures. This is particularly apparent in the c-orbital contribution, which exhibits a significant spectral weight for |k| k F below the Fermi energy at T = 0 [see Fig. 2(a)]. If T gets larger than t c , the c-orbital band leaks into the region above Fermi energy [ Fig. 4(c)] and, as a consequence, c and f bands nearly follow the non-interacting band dispersions α (k) = −2t α cos k with α ∈ {c, f }, see Fig. 4(d). While ARPES on Ta 2 NiSe 5 has not indicated such a behaviour so far, t-ARPES experiments detect the almost bare-band structure showing spectral weight even above Fermi energy [19]. This further supports that the EI state observed in Ta 2 NiSe 5 is of BCS type. Note that the flat valence band shifts to higher energies (but is still located below the Fermi energy) when the temperatures is raised [13]. Since we observe no difference between the location of the flat band near k = 0 at T = 0 and T /t c = 0.1, see Fig. 2(d) and Fig. 4(a), respectively, this might be attributed to the higher dimensionality of the real material.
Let us finally comment on the behavior of the finite-temperature PES in the strongcoupling regime (U/t c = 8; see Fig. 5). Here, by and large, the situation is reminiscent to that in the Hubbard model for U/t = 8, cf. Appendix B.2. While, at the low temperature T /t c = 0.1, the PES barely changes compared to T = 0, at about T /t c = 0.5, a noticeable part of the spectral weight is redistributed to (large) negative energies [see Fig. 5(b)]. As the temperature is further raised up to T J eff = 4|t f |t c /U , a PES band appears above the Fermi energy [ Fig. 5(c)]. It becomes more distinct for T > J eff [ Fig. 5(d)]. Note that in ARPES experiments on Ta 2 NiSe 5 none of these effects have been reported.

Conclusions
To sum up, we numerically examined the ground-state and spectral properties of the extended Falicov-Kimball model (EFKM) in one spatial dimension for the half-filled band case. Using the time-dependent density-matrix renormalization group technique in combination with the purification method and infinite boundary conditions, we derived unbiased results that hold in the thermodynamic limit, i.e. for the infinite system at both zero and finite temperature.
Regarding the ground-state phase diagram we confirmed the existence of an excitonic insulator state, sandwiched between staggered-orbital-ordering and band-insulator phases in the Coulomb-interaction orbital-splitting plane. The excitonic insulator typifies a BCS (Bose-Einstein) macroscopically coherent quantum state for weakly (strongly) coupled c and f electrons. Our analysis of the pair condensation amplitude allows to quantify the BCS-BEC crossover region, where the character of the constituting electron-hole pairs changes from "Cooper-like" (on the semimetal side) to tightly-bound excitons (on the semiconductor side).
The spectral properties of the EFKM are further evidence for a BCS-BEC crossover of the excitonic condensate. Hallmarks of the BCS-BEC crossover are an increasing band renormalization, including a widening of the charge gap, as well as spectral weight transfer (also from the coherent to the incoherent part of the single-particle spectrum) and a weakening of Fermi surface effects.
Most interesting in view of recent ARPES experiments on Ta 2 NiSe 5 seems to be the band flattening, which has been taken as strong indication for the formation of an excitonic insulator state. Simulating the zero-temperature photoemission in the weak-coupling regime of halffilled EFKM with a parameter set that (independently) has been shown being most suitable for Ta 2 NiSe 5 , the flat valence band is confirmed in the narrow region of momentum space detected by ARPES. Therefore we believe that Ta 2 NiSe 5 typifies rather a BCS-type excitonic insulator. Findings that Ta 2 NiSe 5 behaves semiconductor-like in some circumstances might be because it is located close to the band insulator boundary in the EFKM model parameter space.
At finite (high) temperatures, a signature of the conduction band appears in the photoemission spectra above the Fermi energy. Here, both valence and conduction bands follow a rescaled cosine dispersion with different hopping amplitudes. It would be interesting to estimate the effective temperature after pulse irradiation by comparing our numerical data with experimental results. Pulse irradiation can be treated within the t-DMRG scheme as demonstrated in Ref. [45], i.e., we can carry out a real-time evolution to simulate A(k, ω) in this case. Such a direct comparison with time-dependent ARPES experiments on Ta 2 NiSe 5 is left for further studies.

A Numerical approach
In this Appendix, we explain how to simulate the time-dependent correlation functions A ± (x, t) of Eq. (5) by combining the t-DMRG technique with the purification method for matrixproduct states (MPS).
The density matrix ρ(β) of the system can be regarded as the reduced density matrix of a pure state |ψ(β) in an enlarged Hilbert space, ρ(β) = Tr Q |ψ(β) ψ(β)|, where the trace is taken over the space Q spanned by the auxiliary sites. The expectation value of an operatorÔ therefore becomes Ô β = ψ(β)|Ô|ψ(β) . To determine the equilibrium density matrix at target temperature T , we first construct a state |ψ(0) corresponding to ρ(β = 0) and then perform an imaginary time evolution |ψ(β) = e −βĤ/2 |ψ(0) on the physical subsystem. For |ψ(0) in the grand-canonical ensemble, we choose a state with a simple MPS representation, in which each physical site is in a maximally entangled state with an auxiliary site. The time evolution is then carried out, e.g., by using the time-evolving block decimation technique [47,48] and swap gates [49].
To obtain a high momentum resolution for the spectral functions, the correlation functions need to be calculated up to long enough distances. It is therefore important to consider large systems or even work directly in the thermodynamic limit by using infinite MPS (iMPS) with a translationally invariant unit cell. Following the latter approach, we calculate the purification |ψ(β) of the equilibrium density matrix using the infinite time-evolving block decimation technique [47,48] with reorthogonalization [50]. Similar to the zero-temperature case, dynamic properties are calculated by switching to real-time evolution after a local perturbation is applied to |ψ(β) . The perturbation lifts the translation symmetry of the state so that it can no longer be described by an iMPS with a repeated unit cell. It is sufficient, however, to update only the tensors of the MPS in the finite region over which the perturbation spreads during the time evolution. The resulting algorithm is essentially that for a finite-system with specific infinite boundary conditions (IBC) [42].
To extend the simulated time and thereby increase the resolution in energy space, we furthermore exploit the time-translation invariance as described in Refs. [51,52]. The idea is to consider two statesÔ j |ψ(β) andÔ j+x |ψ(β) , and evolve one forward and one backward in time so that their scalar product yields the correlation function for distance x. When doing this, it is important that in addition to the time evolution of the physical sites, the auxiliary sites are evolved in the reverse direction, which also slows down the buildup of entanglement in the purification state [53]. An advantage of IBC, in addition to the absence of finite-size effects, is that we can exploit the spatial translation symmetry of the equilibrium state when calculating the correlation functions with the above method. As demonstrated in Ref. [54], the correlation function at an arbitrary distance can be obtained from a single state just by shifting the two states with forward and backward time evolution relative to each other, while one would need a separate simulation for each distance with open boundary conditions. Special care should be taken in implementing this procedure if the operatorÔ j in Eq. (5) is fermionic as in the calculation of the (I)PES. The Jordan-Wigner stringsF j = exp(iπn j ) withn j =n j,↑ +n j,↓ , which appear in the mapping of spinful fermionic operators into spinful bosonic operators, need to be taken into account when preparing the MPS for the real time evolution. We defineĉ whereâ j ,σ obey fermionic anticommutation relations if j = j and σ = σ , and bosonic commutation relations otherwise. ForÔ j =ĉ j,↑ (Ô j =ĉ j,↓ ), we have to apply the mapped operator j−1 =1F Ô j , whereÔ j =â j,↑ (Ô j =F jâj,↓ ). Note, that the Jordan-Wigner strings only act on the physical sites.
In the following, we give the step by step explanation of the complete algorithm (for simplicity at T = 0): 0) Simulate the ground state |ψ 0 in iMPS representation using iDMRG [55]. Here, we consider a two-site unit cell iMPS in the canonical form with tensors Λ [n] and Γ [n] (n ∈ {0, 1}).
1) Construct an MPS with IBC with appropriate window size L W by repeating the iMPS. The resulting state can be written as where σ j denotes the basis states of the local Hilbert space at site j. After step 4 we go back to step 3 until the target time is reached. At finite temperatures this process can be performed in an analogue manner. Finally, obtained data can be extrapolated to longer times through linear prediction [56]. We apply an exponential windowing function to the obtained time-dependent correlation functions, which, after the Fourier transform, corresponds to a convolution of the spectral functions by a Lorentzian η L /(π(ω 2 − η 2 L )) with the broadening parameter η L . In Appendix B we will exemplarily calculate the spectral functions in the Hubbard model by means of the presented t-DMRG-based scheme. Again the time step τ = τ β = 0.01 and the maximum bond dimension is 800 so that the truncation error is less than 10 −6 .

B Spectral functions in the Hubbard model B.1 Model Hamiltonian
The 1D Hubbard model readŝ whereĉ † j,σ (ĉ j,σ ) creates (annihilates) a fermion with spin projection σ ∈ {↑, ↓} at lattice site j, andn j,σ =ĉ † j,σĉj,σ . Remarkably, at zero temperature, the model (8) can be solved by the Bethe ansatz [57]. Let us again consider the half-filled case where the number of particles N is equal to the number of lattice sites L. In this case, the model has a Mott insulating state for any U > 0 + with a finite charge gap that increases exponentially with U in the weak-coupling regime and grows linearly for large interactions. The spin-degree-of-freedom excitations, however, are gapless, so that the spin-spin correlations decay with a power-law.
The momentum-resolved spectrum of the physical excitations can be obtained by considering two different types of elementary excitations: gapped spinless excitations carrying charge ±e called holons (h) and antiholons (h), respectively, and gapless charge-neutral excitations  carrying spin ±1/2 called spinons (S z = 1/2 spinon s and S z = −1/2 spinons). At half filling, the Bethe ansatz yields the following dressed energies and momenta for these excitations in the thermodynamic limit [57]: where J n (ω) are n-th-order Bessel functions. The dispersion relations E(P) are obtained by varying the parameters Λ ∈ (−∞, ∞) and k ∈ (−π, π] . The physical excitations follow from permitted combinations of the elementary excitations. For example, the spin-charge scattering states whose energy and momentum are given by In the following section B.2, our numerical approach will be benchmarked against the Bethe ansatz results (9)-(11).   Figure 7 presents the PES, A − (k, ω), for a half-filled Hubbard chain in the strong-coupling regime with U/t h = 20 (upper panels) and 8 (lower panels) at T = 0. Data are obtained by the t-DMRG technique with IBC. We observe two separately dispersing peaks in the momentum interval [0, π/2] that merge at the Fermi momentum k F = π/2 where the gap opens. In the interval [π/2, π] we also find two dispersive peaks merging at the zone boundary. Comparing the course of the t-DMRG peak positions with the exact Bethe ansatz dispersions (9)-(11), we are able to identify the physical nature of each dispersion, see right panels of Fig. 7. Since agreement is excellent, not only in the limit of large Coulomb interaction U/t h = 20 but also for U/t h = 8, we can clearly assign the excitations as spinon (red lines) or holon (blue lines) branches, and can also identify the onset of the holon-spinon continuum (green lines). As an advantage, the t-DMRG data also provide the spectral weight as a function of momentum (cf., the color code). Because the spinon bandwidth is about the effective exchange coupling J eff = 4t 2 h /U , it is almost flat for very large U and gets wider when U shrinks (cf. upper and lower panels). We would like to point out that the upper onset of the peaks of A − (k, ω) nicely agree with half of the single-particle gap ∆ c /2 , see white and black lines. The whole single-particle spectrum A(k, ω), which includes also the IPES contribution, is shown Fig. 8 for interaction strengths ranging from U/t h = 1 to 8. Since the Mott gap ∆ c is exponentially small in the weak-coupling regime [∆ c /t h 0.005 (0.173) for U/t h = 1 (2)], the dispersion for U/t h = 1 [ Fig. 8(a)] is practically identical to that of the bare cosine band with the holon bandwidth W h = 4, without opening a gap due to the finite Lorentzian broadening (η L /t h = 0.1). Increasing U slightly, both spinon and antiholon branches appear, see Fig. 8 . We note that A(k, ω) for U/t h = 4 was previously discussed within t-DMRG in Ref. [44], but only for finite systems. Our infinite-system data validate these results to a large extent. Figure 9 demonstrates the temperature effects on PES in the half-filled 1D Hubbard model at various Coulomb interaction strengths obtained by combining the t-DMRG technique with the purification method as described in Appendix A. Clearly, at low temperature (T /t h = 0.1) the spectral functions are very similar to their zero-temperature counterparts, whereby the peak maxima are slightly reduced (see the different scales of the colorbars with increasing temperature). In the strong-coupling regime, the two-brunch structure of the spectral function only survives up to T ∼ J eff (cf. the panels for U/t h = 4 and 8 in Fig. 9). Furthermore, the "V-shaped" structure of the spectrum and the concentration of spectral weight around ω = −U/2 − 2t h , k = 0 and ω = −U/2 + 2t h , k = ±π found in the strong-coupling approach [44] is corroborated by our numerical results. A signature of upper Hubbard band is also visible here, which can be qualitatively described by the Hubbard-I mean-field approximation [58,59] (note that we have shown only the PES part of A(k, ω) in Fig. 9).

B.3 Photoemission spectra at finite temperatures
Most interesting seems to be the melting of the Mott gap in the vicinity of the Fermi level: The zero-temperature PES spectral weight leaks into the Mott gap regime in the doped Mott insulator (two holes in the finite-size system), which leads to the extension of the spinonantiholon continuum [44]. Similarly, holes in the lower Hubbard band left behind as a result of doublon creations can be occupied by thermally excited quasiparticles. This leads to the melting of the Mott gap with increasing temperature, which becomes more distinct when decreasing the Coulomb repulsions U .
At very high temperatures (T > t h ), the energy distance between upper and lower Hubbard bands is insignificant (provided U is not too large) and finally both bands merge into a single band (k) = −2t h cos k, see Fig. 9(p). k −π 0 π k −π 0 π k −π 0 π k −π 0 π Figure 9: Photoemission spectra A − (k, ω) in the half-filled 1D Hubbard model for various Coulomb repulsions U/t h at different temperatures T /t h .