In-Gap Band Formation in a Periodically Driven Charge Density Wave Insulator

Periodically driven quantum many-body systems host unconventional behavior not realized at equilibrium. Here we investigate such a setup for strongly interacting spinless fermions on a chain, which at zero temperature and strong interactions form a charge density wave insulator. Using unbiased numerical matrix product state methods for time-dependent spectral functions, we find that driving of the correlated charge-density wave insulator leads not only to a renormalization of the excitation spectrum as predicted by an effective Floquet Hamiltonian, but also to a cosine-like in-gap feature. This is not obtained for a charge density wave model without interactions. A mean-field treatment provides a partial explanation in terms of doublon excitations. However, the full picture needs to take into account strong correlation effects.


I. INTRODUCTION
A central driving force of modern condensed matter physics is the realization of novel phases of matter out-of equilibrium like transient superconductivity [1][2][3][4][5] or excitonic insulators in transition metal dichalcogenides [6][7][8][9][10][11][12]. These are typically created from a highly complex interplay of band structure, (electronic) interactions, and excitations by a light-field. Nowadays, experimental techniques allow to actively "engineer" properties of many-body quantum systems in such out-of-equilibrium systems in a highly controlled way [13][14][15]. This opens up the possibility to create behavior that is not even possible in equilibrium setups. An important pathway is to induce excitations whose interplay with the electronic interactions can lead to intriguing transient behavior. Such excitations can be realized in experiments, e.g., by ultrashort laser pulses in so-called pump-probe setups [16][17][18], or by continuous periodic driving of the systems, e.g., by shaking ultracold atom systems on optical lattices [19][20][21][22][23][24][25]. A much studied theoretical idealization is (infinite) periodic driving, which can be addressed by Floquet theory. In this framework the properties of the system can be described by a time-independent effective Hamiltonian [26,27]. Control over the parameters of the driving translates to control over the effective Hamiltonian and this, in turn, allows to manipulate order parameters [28], induce topological order [29] etc.
It is possible to derive time-independent approximate effective Hamiltonians in the low- [30] and highfrequency [26,31] regimes, and development of methods in this vein is ongoing [32]. In addition, the periodic driving generically leads to energy absorption evolving the system towards an infinite-temperature state [33,34]. This restricts the relevance of such effective Hamiltoni-ans to regimes in which the energy absorption is suppressed like the high-frequency regime. Additional transient dynamics can be induced by the switching-on procedure [27,[35][36][37][38]. One major topic of interest is the role of interactions in strongly driven systems. It has been shown that in Hubbard systems at resonance the interaction can be renormalized, and double occupancies can be enabled [26,31,39]; one can even tune the parameters of the driving so that the fermions behave like free particles [40].
In our work, we investigate driven strongly interacting fermions including a sudden switching-on procedure and without assuming the high-frequency approximation. This allows us to look for new effects beyond these limiting cases. Therefore, we calculate non-equilibrium spectral functions with unbiased matrix product state (MPS) [41,42] approaches, and systematically investigate correlation effects by comparing to non-interacting and mean-field scenarios. This is in contrast to other approaches, which e.g. rely on the Floquet-Magnus expansion [26,40] in terms of the inverse driving frequency. We focus on a simple paradigmatic model for strongly correlated physics, namely a chain of spinless fermions with nearest-neighbor interactions. At half filling and zero temperature the model is known to undergo a Berezhinskii-Kosterlitz-Thouless (BKT) type transition [43] from a Luttinger liquid (LL) [44] to a interacting charge density wave (CDW) insulator [45] when increasing the interaction strength. Driving the system with frequencies much larger than the gap ("Magnus case"), a renormalization of the gap size for this system is predicted [28]. We investigate this by calculating non-equilibrium single-particle spectral functions. For a static band structure upon monochromatic periodic driving, one expects on general grounds the formation of Floquet side-bands [46]. Our approach allows us to study both effects, but also to go beyond and to look for the emergence of additional spectral features.

II. DRIVING A STRONGLY CORRELATED CHARGE-DENSITY-WAVE INSULATOR
We consider a periodically driven chain of interacting spinless fermions described by the Hamiltonian where V is the strength of the density-density interaction and µ i is an on-site potential. A(t) = θ(t)A 0 sin Ωt is a time-dependent vector potential, which is switched on at time t = 0. This type of coupling is known as Peierls substitution [47] and models a monochromatic classical light-field, which couples to the electrons in the system (here: spinless fermions). We will mostly consider open boundary conditions (OBC) and for comparison periodic boundary conditions (PBC). For A(t = 0) Bethe ansatz (BA) [48] gives the BKT transition at V /t h = 2. In the following, we will consider driving a system in the CDW phase at V /t h = 5, for which the energy gap according to BA [48] is ∆/t h ≈ 1.576. We define the non-equilibrium generalization of the spectral function via the Fourier transform of the retarded Green's function with a damping factor η ≈ 0.1 as further explained in App. B. Integration of A ret k (t, ω) over the crystal momenta k directly yields the time-dependent density of states (tDOS). This quantity relates to measurements in time-dependent angle resolved photoemission spectroscopy (trARPES), although a more detailed modelling is required for a direct comparison to experiments [17,49,50], as well as for its interpretation in the deep nonequilibrium regime [51]. Nevertheless, these quantities allow us to study qualitative changes of the spectral function with time, such as the formation of additional branches, or a change of the band structure due to the excitation.
We complement our study by considering the time evolution of the CDW order parameter further details are explained in App. B. We keep also track of the time evolution of the energy, E(t) = H(t) , which serves as a measure for the heating of the system.

III. RESULTS
For periodic driving well above the band gap ("Magnus case") the properties of the driven model are expected to be well-described by an effective Hamiltonian according to the Floquet-Magnus expansion [28,31]. It is given by the original Hamiltonian but with a renormalized hopping parameter t eff h = J 0 (A 0 )t h ≈ 0.7652t h for A 0 = 1 as in Ref. 28. This corresponds to a model with V /t eff h 1.3V /t h . In the following, we consider these aspects by studying the system at V /t h = 5, which is deep in the CDW insulating phase.
Figs. 1 and 2 show results for the equilibrium and for the non-equilibrium spectral functions at waiting times t = 0, 10t −1 h , and 20t −1 h . The driving frequency of Ω/t h = 10 > ∆/t h ≈ 1.576 is larger than the spectral gap but not substantially larger than the interaction strength so that we are not deeply in the Magnus regime. Fig. 1(e) shows the equilibrium spectral function for the effective Hamiltonian with the renormalized hopping matrix element. Let us consider the equilibrium results first: In Fig. 1(a) we can identify the equilibrium continuum of excitations [52] as well as the spectral gap located around ω = 0. Despite finite size effects its minimal size around k = π/2 is in agreement with the BA prediction. The spectral function of the effective Hamiltonian in Fig. 1(e) looks similar to Fig. 1(a) but the width of the continuum is smaller and the gap is larger. Turning to the nonequilibrium results, we see that the spectral function at waiting time t = 0 in Fig. 1(b) looks quite similar to the equilibrium result for the effective Hamiltonian in Fig. 1(e), but it possesses additional features. In particular, a new in-gap band comes into appearance and the continuum changes slightly its size and form. Around frequencies ω/t h = ±10 weak signals appear, which seem to echo the in-gap feature. These are reminiscent of Floquet sidebands, which are observed in time-resolved ARPES experiments [46]. It is noteworthy, however, that one seems to obtain these "echoes" only for the in-gap signal, but not for the main spectral features. At later waiting times t = 10t −1 h and t = 20t −1 h , the shape of the continuum does not further change, but the in-gap signal becomes more pronounced. This is further confirmed by looking at Fig. 2, which shows a momentum cut through the central region of the spectral function at k ≈ π/2 for waiting times t = 0, 10t −1 h and 20t −1 h . In order to better focus on the relative distribution of spectral weight we normalized all spectral functions to their maximum value at that k-slice. In all cases, we identify two main lobes and smaller peaks. Let us follow the behavior of the main lobes and of the largest peaks at ω = ±10 and ω = 0: At waiting time t = 0, the lobes show a small difference between the nonequilibrium result and the result of the effective Hamiltonian, which becomes even smaller at later waiting times. It is noteworthy that the similarity to the effective description is already obtained at waiting time t = 0, although the effective Hamiltonian relied on the infinite-driving assumption. The additional h . The equilibrium spectral features in e) are also present in the spectral function of the driven system. In addition, there is clearly additional spectral weight appearing whose main feature is well approximated by f (k) ≈ −2.5 cos(ka) (dashed line in d)). All data is obtained with MPS time-evolution for a system with L = 64 chain sites and open boundary conditions.  Cross section plots at k ≈ π 2 through the data in Fig. 1. The dotted grey lines present the equilibrium spectral function of the undriven system, the dashed black lines correspond to the system with Floquet-renormalized parameters. The data shown is normalized to the maximum value at this k-slice:Ã k (t, ω) = A k (t, ω) / maxω A k (t, ω). The solid lines correspond to the non-equilibrium spectral functions (red color = positive, blue color = negative). signals at ω/t h = ±10 oscillate in time, but are suppressed with increasing time. Nevertheless, on the time scale treated by us, the peak at ω = 0 becomes more pronounced with time and so the cosine-like in-gap feature of Fig. 1 appears stable on the time scales treated by us. We checked that it is also present for other system sizes and periodic boundary conditions (PBC) so that boundary effects can be ruled out as an explanation (cf. App. A). By comparing results for L = 32 and L = 64 we find that the peak gets sharper for larger system size, while keeping the relative weight. At early waiting times negative weight appears in the spectral function. This is not an artefact and traces back to the non-equilibrium nature of the state. It was reported recently [36,53] that upon averaging of the Wigner average time coordinate over a driving period, the non-equilibrium density of states for fermions can be shown to be positive. However, at later waiting times away from the turning-on of the field at time t = 0, our spectral function (obtained using "horizontal time coordinates" [36]) is also almost completely positive without this procedure.
To better understand our findings, we study in Fig. 3 the system at the same value of V /t h = 5 but with a driving frequency closer to resonance Ω/t h ≈ 42 · 10 −1 for comparison. The additional feature in the gap region in this case is even stronger pronounced and it goes hand in hand with a significant reduction of the original spectral features of the CDW insulator. The question arises how this disappearance of the quasiparticle continuum is connected with a destruction of the CDW state. To study this, we calculate the time-evolution of the CDW order parameter O CDW (t), which is displayed in Fig. 4 for driving in the Magnus regime and closer to resonance. In the latter case, O CDW (t) completely vanishes on a time scale t ≈ 10t −1 h , which is in agreement with the time scale on which the holon continuum disappears in the spectral function. The behavior for Ω/t h = 10 is more complicated, but also here a partial melting of the CDW state is realized, which continues over times longer than the ones treated by us. Note that the parameters of the effective Hamiltonian are deeper in the CDW phase and in equilibrium one would expect a larger CDW order parameter. In contrast, here we observe a melting of the order, which is due to the nonequilibrium protocol applied and the absorption of energy. Driven systems in the long-time limit will realize an infinite-temperature state [33,34]. In our case, as can be seen in Fig. 4(b), the energy continues to increase as a function of time indicating that, on the transient time scale treated by us, the infinite temperature state is not yet realized. Clearly, energy absorption is increased closer to resonance.
We would like to distinguish our finding further from a known effect: In earlier works on electron-mediated CDW melting [54,55], the appearance of in-gap spectral weight was reported already in a pumped noninteracting fermion model as a genuine non-equilibrium effect. Hence, the question arises, if the new in-gap band can be obtained also in a continuously driven CDW system without interactions. To study this, we adopt the "A-B model" by Shen et al. [54] at half filling, and apply the same semi-infinite driving protocol used for the t-V chain. The CDW order in the model is due to the presence of a staggered on-site potential and leads to a spectral gap of size ∆ ≈ U . We use the same Trotterized time-evolution as in the original work [54] and choose a step size of 10 −6 t −1 h . The results of the simulations for a driving frequency of Ω = 10t h and a gap of U = 5t h are shown in Fig. 5. The dynamics of the order parameter in Fig. 5(d) is more oscillatory than in Fig. 4, and its envelope is decreasing in time, indicating CDW melting. The momentum cuts through the spectral function, however, show that this is not connected with the formation of a peak in the spectral gap. One should note that the spectral function does not become stationary in the model but in the in-gap region the only effect appears to be a small shift (compare Fig. 5(g) and (h) ).
To go beyond the purely non-interacting limit we treat the dynamics of the CDW phase in the t-V model within a Hartree-Fock mean-field (MF) approach, whose results are shown in Fig. 6. This will allow us to investigate the role of interaction-induced doublon excitations for the in-gap feature. A more detailed discussion can be found in App. C. The equilibrium MF band structure is similar to the one in the A-B model. In the driven model, however, we obtain, in addition to the Floquet replicas of the equilibrium bands, a signal in the band gap around ω/t h ≈ 0. It stems from additional resonances of the band structure at separation ∆ω = ±V , which come from the breaking of a doublon present in the system, or the formation of a doublon, respectively (a doublon here is formed by two electrons on adjacent sites). The larger V , the fewer doublons are present in the ground state, so that the resonance at −V is weaker, as confirmed within our MF approach (see App. C).
This picture will explain in parts the MPS findings.  However, they differ from the MF results in at least two aspects: i) the MF results cannot reproduce the exact quasi-particle continuum and hence the in-gap resonance possesses further features caused by higher scattering processes or correlation effects. Also, more complicated excitations can come into play, e.g., bound states, which are not captured by MF. ii) when going closer to resonant driving, the melting of the CDW rather seems to induce a new band than a replica of existing spectral features, see Fig. 3. Taking these considerations into account, we propose the following possible scenarios leading to the in-gap band observed by MPS: i) creation of bound states described in Ref. 52 for the system at equilibrium. This would lead to a cos-like band whose bandwidth, however, depends on the interaction strength ∼ 1/V . The findings of Figs. 1 and 10 for V /t h = 5 and V /t h = 10, respectively, show that the bandwidth in our case seems to be only weakly depending on V , if at all.
ii) Formation of doublons and scattering of the elementary excitations. A reminiscent scenario is, e.g., realized in spin-chains upon increasing the temperature [56,57].
To investigate this, one should treat the single-particle excitations of the system, e.g., with Bethe ansatz, which is beyond the scope of this paper.
iii) the melting of the CDW state leads to the 'emission' or creation of free carriers, which can move freely on the lattice and hence realize a tight-binding like dispersion. In this scenario, the incomplete melting of the CDW would lead to a situation in which remnants of the CDW crystal survive, but at the same time the system would possess also some metallic character due to the freely mobile charge carriers.
The comparison with the MF-results indicates that scenario ii) is probably best applicable deeper in the Magnus regime, while closer to resonance scenario iii) might be better suited. Typically, a mix will be realized. In addition, as indicated by Figs. 7, 8 and 10 for various parameters, there seems to be a broadening or a continuum attached to the cosine-like feature. Further investigations are needed to clarify this.

IV. CONCLUSIONS AND OUTLOOK
We observe by computation of the time-dependent spectral function using MPS that the melting of an initial CDW insulator is accompanied by the formation of a new cosine-like feature in the spectral function upon periodic driving with frequencies above the band gap ("Magnus regime"). When approaching resonance, the non-equilibrium spectral function changes significantly as compared to the equilibrium case and the new band becomes the dominant feature. Instead, at frequencies substantially larger than the gap the original features of the spectral function are modified according to the prediction by the effective Floquet Hamiltonian [28], but in addition the cosine-like in-gap band prevails. Such a feature is not observed when periodically driving a noninteracting CDW-state. On the MF level, an in-gap band is obtained by creating or breaking a doublon, but its properties still differ significantly from the MPS results, indicating that correlation effects are needed to explain our findings. It will be interesting to further clarify the origin of this new band in this simple correlated system, e.g., using higher-order Floquet-Magnus expansions, Bethe ansatz, or semiclassical approaches, such as fermionic truncated Wigner approximations [58].
It is an open question to see whether the interplay of the melting of CDW states and electron correlations can lead to similar in-gap features beyond the leading order Floquet-Magnus effective Hamiltonian also in periodically driven interacting twodimensional systems, such as tilted bilayer heterostructures.

ACKNOWLEDGMENTS
We thank Götz Uhrig, André Eckardt, Stefan Kehrein, Sebastian Paeckel, Thomas Köhler, Mona Kalthoff, Karlo Penc, Niklas Bölter and Karun Gadge for useful discussions. We are grateful for many stimulating and insightful discussions with all participants of the journal club of the B07 project of the SFB 1073, in particular also Stefan Mathias and Marcel Reutzel. The work was supported by the North-German Supercomputing Alliance (HLRN) and we are grateful to the HLRN supercomputer staff. We also acknowledge access to computational resources provided by the GWDG and acknowledge technical assistance. This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -217133147/SFB 1073, projects B03 and B07. The MPS-based results presented in this work were generated using the SymMPS toolkit [59].

Appendix A: Additional supporting figures
In this Appendix we present some additional figures with other parameters and more technical content. Figs. 7 and 8 provide additional cross section plots at momenta k ≈ 0 and k ≈ π/4. One can clearly identify the in-gap spectral weight as well. It is noteworthy that here it is less sharply concentrated and a continuum of weight seems to be attached to the peak.
In Fig. 9 we collect data for different system sizes and varying boundary conditions as a check that our observations do not depend on these parameters. Looking at Fig. 9(g) to (i) the peak at frequency zero is clearly visible for all data sets. However, for L = 32 and OBC it is less pronounced. The PBC data and OBC with 64 sites agree very well with each other. Fig. 10 shows a simulation with a larger value of V = 10t h and double the driving frequency Ω = 20t h . Our observation from the main text is confirmed and  we find the cosine-like feature again. It follows roughly the same functional form as the signal in the main text although here the interaction strength is different (see Figure caption) Appendix B: Details on Green's functions, spectral functions

Definitions
All Green's functions (GFs) are derived from the contour-ordered single-particle Green's function [60] which can be written in a matrix representation with respect to the forward and backward branches of the realtime axis. In this representation the greater G > αβ (t, t ) and lesser Green's function G < αβ (t, t ) each have one time argument lying on the forward and one on the backward branch of the real-time contour. The retarded GF is a linear combination of the two with an additional theta function,  At equilibrium, one of the two time variables can be suppressed due to time-translational invariance. Outof-equilibrium, however, we need to consider both time variables, and the Fourier transform to frequency space is not unique any more (see, e.g., Ref. 36). In order to minimize the numerical costs, we choose to use "horizontal" time coordinates, in which we evolve the wavefunction up to a time t and then perform the Fourier transform with respect to the relative time τ = t − t (also referred to as waiting time farther below), after further evolving the system in time, with t > t The states are labelled by momentum indices k (depending on boundary conditions, see below). In the numerics we calculate the auxiliary quantities such that and we finally obtain the nonequilibrium spectral functions

Interpretation
In equilibrium the retarded GF contains information about the density of states while the lesser and greater GFs contain information about occupations of the states. The latter is reflected in the fact that due to the absence of the θ(τ ) the whole "history" of the lesser and greater GFs needs to be traced back until time −∞. In this text we compute the retarded Green's function. Integration of A ret k (t, ω) over k directly yields the time-dependent local density of states. In contrast to the equilibrium sitatuation the quantity A ret k (t, ω) is not necessarily positive, so that care needs to be taken with this interpretation [51]. Further issues when using the two-time Green's functions to directly model trARPES experiments arise due to the lack of gauge invariance, in particular if pump and probe pulses overlap, or if the system possesses multiple bands, see Refs. 17, 49, and 50 for a detailed discussion. Here, we neglect these aspects and mainly focus on the most important qualitative features evolving with time. However, central predictions of Floquet theory, like the effective Hamiltonian picture, are captured with a good accuracy, so that we believe that our results not only provide a qualitative picture, but also quantitative predictions with a good accuracy.
Temporal Fourier transform Due to the finite maximal τ that we are able to reach and due to the theta function, the temporal Fourier transform produces a nonzero background ∼ 10 −5 signal everywhere in the spectral function. Since we expect the exact spectral function to have value zero if no signal is present, we decided to subtract this background. Note that we add a damping factor η ≈ 0.1 to regularize the finite time-propagation with respect to τ . In addition, to improve the ω-resolution we padded the τ -data with zeros to obtain at least 4096 frequency points.
Symmetry-broken CDW state For MPS calculations with OBCs we study the system at half filling with µ i = 0. In the CDW phase, this leads to an exact superposition of the two possible symmetry-broken ground states, so that for a finite system O CDW (t) = 0 for all times t. In order to be able to keep track of the dynamics of the CDW order parameter, we performed additional simulations where we applied a 'pinning field' at the edge µ 1,L = 0, which selects one of the ground states and allows us to study O CDW (t). We have checked that the spectral function does not differ much if calculated with  Fig. 1 The data is obtained for OBC in a chain of L = 32 sites.
or without the pinning field. In order to minimize possible boundary effects in the order parameter, we perform the sum in Eq. (3) only over the four unit cells in the center of the system.
MPS calculations In order to compute G ret,≶ kk (t, τ ) using matrix product states we start from the system's ground state |GS and consider the following quantum states where U TDVP (t, t 0 ) denotes time-evolution using the MPS realization of the time-dependent variational principle (TDVP). Here, we always apply a two-site TDVP algorithm [42]. The operators c k carry momentum space labels. However, we always work in position space and exploit that we can write momentum space annihilation and creation operators as a sum of local operators c k = j P k,j c j , where we have introduced the transformation matrix P , allowing us to compute C ≶ kl (t, τ ) through a series of local MPO-MPS applications. P depends on the boundary conditions used. Using the states (B6) we calculate the quantities which are related to G ret,≶ kk (t, τ ) as desctibed in the Definitions section. For PBC we implement a "snake geometry" [63] for the labelling of the sites in the chain. It turned out that the DMRG ground state search always chose one of the degenerate ground states which allowed us to calculate the order parameter directly. Further de-tails on the calculation of the non-equilibrium spectral function with MPS methods can be found in Ref. 62.

Appendix C: Hartree-Fock time-evolution
We start from a Hartree-Fock decoupling of the interaction term and assume a two-site unit cell with sublat-tices A and B. Let us denote ρ A := c † i c i i∈A , ρ 0 := c † i c i+1 i∈A , ρ B := c † i c i i∈B , ρ 1 := c † i c i+1 i∈B . (C1) Using the Fourier basis (Q = π) we obtain, using the definitions k = −2t h cos(k), χ k = V ρ 0 e −ik + ρ * 1 e ik , the following representation of the Hamiltonian In the following we consider half filling µ = V (ρ A +ρ B ). The saddle point values of ρ A , ρ B , etc. are determined with a simulated annealing approach. Diagonalization of the Hamiltonian yields the eigenenergies Hence, the spectral gap is given by 2V (ρ B − ρ A ). For the dynamics we first solve the time-diagonal problem and obtain the full one-particle reduced density matrix ρ ij (t) = c † i (t)c j (t) . In a second iteration we solve the equation of motion for the relative time τ using the time-diagonal data from the first iteration. This corresponds to solving the Kadanoff-Baym equations with a Hartree-Fock self-energy [60].
We have studied chains with L = 64 sites and periodic boundary conditions for which we calculated the timeevolution of the driven and undriven model. For V = 8t h (chosen to have a visual separation of V -dependent spectral features from Floquet sidebands) the MF spectral gap size is about ∆ ≈ 15t h . We choose a driving frequency of Ω = 30t h , which yields the same ratio ∆/Ω as for the A-B model. The results are shown in Fig.  6. The undriven spectral function is very similar to the one obtained in the A-B model although with a different spectral gap size. The driven model, however, shows additional signals separated by ±V from the main peaks. This gives rise to an in-gap signal around ω/t h ≈ 0. Like in the A-B model the spectral function displays negative weights, which -in contrast to our MPS results -FIG. 11. Simulation results for the t-V model within Hartree-Fock mean-field theory with V = 12t h and for semi-infinite sinusoidal driving with amplitude one and frequency Ω = 30t h (cf. Fig. 6) The data is obtained for a model with L = 64 sites and periodic boundary conditions. are pronounced at all times treated by us, while in the MPS case the negative weights seem to substantially decrease in time. We consider, however, that it still contains relevant qualitative information, e.g. the position of spectral peaks. The order parameter in Fig. 6(e) is oscillatory with a slowly decaying envelope. It in fact oscillates around a larger order parameter than in equilibrium. Still, the equilibrium CDW is broken up and charges can move in the system. We have checked that when the driving is suddenly turned off order parameter oscillations as well as the in-gap spectral features remain. The Floquet replicas of the main bands, however, disappear.