Exploring nonequilibrium phases of photo-doped Mott insulators with generalized Gibbs ensembles

Many experiments show that strong excitations of correlated quantum materials can cause non-thermal phases without equilibrium analogues. Understanding the origin and properties of these nonequilibrium states has been challenging due to the limitations of theoretical methods for nonequilibrium strongly correlated systems. In this work, we introduce a generalized Gibbs ensemble description that enables a systematic analysis of the long-time behavior of photo-doped states in Mott insulators based on equilibrium methods. We demonstrate the power of the method by mapping out the nonequilibrium phase diagram of the one-dimensional extended Hubbard model, which features η-pairing and charge density wave phases in a wide photo-doping range. We furthermore clarify that the peculiar kinematics of photo-doped carriers, and the interaction between them, play an essential role in the formation of these non-thermal phases. Our results establish a new path for the systematic analysis of nonequilibrium strongly correlated systems. Light-matter interactions can induce phases and physics not observable in the equilibrium regime of a quantum material but understanding the underlying mechanisms can be challenging. Here, the authors show a generalised Gibbs ensemble applied to photo-doped states in Mott insulators can be used to investigate the non-thermal phases of various correlated states.

N onequilibrium control of quantum materials is an intriguing prospect with potentially important technological applications [1][2][3][4][5][6] . Experiments with various materials and excitation conditions have reported phenomena not observable in equilibrium including nonthermal ordered phases such as superconducting (SC)-like phases 4,[7][8][9] , charge density waves (CDW) [10][11][12] and excitonic condensation 13 . Among various nonequilibrium protocols, photo-doping is a basic and important one, in which a radiation pulse creates electron-and hole-like charge carriers with a long lifetime on the electronic timescale. Due to the nonthermal distribution of the carriers and the cooperative interplay between them, novel nonequilibrium phases can be induced.
The theory of photo-doping has been extensively discussed for semiconductors, where a rich phase diagram including electronhole plasmas and exciton gases [14][15][16][17][18][19] as well as exciton condensation [20][21][22] is found. The physical picture is that, after electrons and holes are created, they rapidly relax within the conduction and valence bands, while their recombination occurs on a much longer timescale. Thus, at the single-particle level, the numbers of electrons and holes are separately conserved, so that in the intermediate time regime one has a pseudoequilibrium state that can be described by an effective equilibrium theory with separate chemical potentials for the electrons and holes [17][18][19] .
A situation of great current interest is the photo-doping of Mott insulators. In these systems, exotic equilibrium states such as unconventional SC phases emerge upon chemical doping 23 , while photo-doping creates novel pseudoparticle excitations not easily represented in a single-particle picture, e.g., doublons and holons in the single-band case. When the Mott gap is large enough, these excitations are long-lived due to the lack of efficient recombination channels [24][25][26][27][28][29] . Therefore, as in semiconductors, a fast intraband relaxation results in a long-lived quasi-steady state. Previous studies based on short-time simulations indicated the emergence of enhanced CDW 30 or SC correlations [31][32][33][34][35] , as well as novel spin-orbital orders 36 . However, unlike in the semiconductor cases, the long-time behavior of photo-doped Mott insulators is not well understood, due to the lack of appropriate theoretical frameworks.
Steady-state formalisms in which explicit heat/particle baths or other dissipative mechanisms are attached to the system have been recently applied [37][38][39] . These formalisms, however, require attention to the influence of the baths and dissipations, and the use of explicitly nonequilibrium methods. An alternative approach is a pseudoequilibrium description as in conventional semiconductors. Crucial differences from semiconductors are that the approximately conserved entities are pseudoparticles (local many-body states) not explicitly appearing in the initial Hamiltonian and that the approximate conservation of their number is not manifest in the Hamiltonian but arises from kinematic constraints. Previous work 37,[40][41][42][43] introduced the idea of using the Schrieffer-Wolff (SW) transformation to reformulate the problem in a way that explicitly references the approximately conserved quantities and isolates the terms that eventually lead to full equilibration. However, applications of such effective descriptions have been so far limited to small clusters and weak 40,41,43 or extreme excitation conditions 42 .
In this work, we introduce a generalized Gibbs ensemble (GGE) type description for the effective model obtained from the SW transformation by incorporating different "chemical potentials for pseudoparticles", i.e., for the many-body local states. This effective equilibrium description allows us to systematically scan the nonequilibrium states in photo-doped Mott insulators for extended systems using established equilibrium methods. We use this approach to identify and study emerging phases in the photodoped one-dimensional extended Hubbard model. We determine the nonequilibrium phase diagram, where η-pairing 33,35,37,42,44 and CDW phases appear in a wide doping range and reveal the corresponding spectral features. The CDW phase is strongly favored in photo-doped systems, compared to the chemicallydoped ones, and it is characterized by unbound doublons and holons in contrast to photo-doped semiconductors where electron-hole binding is an important effect. We show that the kinematics of photo-doped doublons/holons, which is qualitatively different from the electron/hole dynamics in conventional semiconductors, plays an essential role and leads to the development of CDW and η-pairing correlations described by squeezed systems without singly occupied sites.

Results
GGE description for photo-doped Mott insulators. The generic formulation of the GGE description is given in Supplementary Note 1. Here we explain the procedure focussing on the extended Hubbard model, whose Hamiltonian iŝ withĤ U ¼ U∑ i ðn i" À 1 2 Þðn i# À 1 2 Þ the on-site andĤ V ¼ V∑ hi;ji ðn i À 1Þðn j À 1Þ the nearest-neighbor interaction.ĉ y iσ is the creation operator of a fermion with spin σ at site i,n iσ ¼ĉ y iσĉiσ the spin-density at site i,n i ¼n i" þn i# , and 〈i, j〉 denotes pairs of nearest-neighbor sites. t hop is the hopping parameter. For large U, the half-filled equilibrium system is Mott insulating.
Photo-doping the Mott insulator creates doublons (doubly occupied states) and holons (empty states). When the Mott gap is large, the number of these excited local states is approximately conserved for kinematic reasons [24][25][26][27][28][29] . However, the original Hamiltonian explicitly contains recombination terms, which generate virtual processes that affect the physics even when recombination is kinematically suppressed. In order to explicitly remove the recombination terms, while effectively taking account of the effects of virtual recombination processes, we apply the SW transformation 37,45,46 U;shift describes the shift of the local interaction.Ĥ 3Àsite represents three-site terms such as correlated doublon hoppings, see Method.Ĥ dh;ex sets the correlations between the neighboring doublons and holons aŝ H spin;ex sets spin correlations between singlons (singly occupied states). In this sense, the above model is a natural extension of the t-J model 23 , which is obtained by assuming that either holons or doublons are added (chemical doping) and ignoringĤ 3Àsite .
Due to intraband scattering and environmental coupling (e.g., phonons), intraband relaxation occurs and the system reaches a steady state. Since the numbers of doublons and holons are conserved in the effective model, the steady state can be described by introducing separate "chemical potentials" for them. The corresponding number operators areN holon ¼ ∑ ini;h withn i;h ¼ ð1 Àn i" Þð1 Àn i# Þ andN doub ¼ ∑ ini;d withn i;d ¼n i"ni# , respectively. With these, the grand-canonical effective HamiltonianK eff can be written asĤ eff À μ holonN holon À μ doubNdoub , or where μ U = μ doub + μ holon and μ = −μ holon . Thus, the local interaction is modified from U by the photo-doping (U − μ U ), i.e., the energy difference between the doublons and holons is effectively reduced, analogously to the effective shift of the band splitting in photo-doped semiconductors 22,47 . The properties of the nonequilibrium steady states may then be described by the density matrixρ eff ¼ expðÀβ effKeff Þ with an effective temperature T eff = 1/β eff 17,22,48 , which is a sort of GGE 49,50 . This is essentially an equilibrium problem that can be studied with established equilibrium techniques. Response functions Àih½ÂðtÞ;Bð0Þ ± i of the nonequilibrium states can also be computed within this framework, see Supplementary Note 2. Our basic assumption that the nonequilibrium states can be characterized by a few parameters, such as the doublon number and effective temperature, is supported by a recent study 37 which demonstrated a good agreement between time-evolving states and nonequilibrium steady states weakly coupled to thermal baths.
Nonequilibrium phase diagram of the photo-doped extended Hubbard model. We apply the above framework to the half-filled one-dimensional extended Hubbard model using infinite timeevolving block decimation (iTEBD) 51 and exact diagonalization (ED) 52 . As in the case of the t-J model, the effect ofĤ 3Àsite is not essential. We confirm this for the photo-doped situation using ED in Supplementary Note 5. Thus, in the following, we focus on the effective modelĤ eff2 , which ignoresĤ 3Àsite . We consider cold systems (T eff = 0) to clarify the possible emergence of nonequilibrium ordered phases. Such a situation may be achieved by energy dissipation to the environment 47,53,54 or entropy reshuffling 55,56 . In the following, we use t hop as the energy unit, and fix U = 10, i.e., the exchange energy is J ex 4t 2 hop U = 0.4. We examine the charge correlations χ c ðrÞ 1 N ∑ i hðn iþr À n av Þðn i À n av Þi, spin correlations η-pairing is characterized by staggered SC correlations. Note thatĤ,Ĥ eff , andĤ eff2 are SU c (2) symmetric with respect to the η-operators for V = 0 44,57 . Due to this symmetry, a homogeneous state with long-range η-SC correlations is on the verge of phase separation, which we avoid by considering nonzero V 42 , see Supplementary Note 4.
In Fig. 1, we show the computed nonequilibrium phase diagram for the photo-doped Mott insulator in the plane of the doublon density (n d ¼ 1 In one-dimensional quantum systems, spatial equal-time correlations can show quasilong-range order, i.e., power-law decay with a critical exponent a less than 2, which corresponds to a diverging susceptibility in the low-frequency limit 58 . The corresponding spatial dependence of the correlation functions is shown on a normal scale in Fig. 2a and on a log scale in Fig. 2b-d. We see that generically more than one correlation function exhibits quasi-long-ranged order. The phase shown in Fig. 1 is identified from the correlation function with the smallest critical exponent. Without photo-doping, an SDW phase with staggered spin correlations is found. However, other correlations quickly become dominant with photo-doping. When V ≲ 0.2 ð¼ J ex 2 Þ, the η-SC phase emerges in a wide photodoping range. This is consistent with recent dynamical mean-field theory (DMFT) analyses for the pure Hubbard model in infinite spatial dimensions employing entropy cooling or heat baths 37,55 . Importantly, the sign of the SC correlations remains staggered regardless of doping and V. For larger V, the CDW phase is stabilized. We note that, in the extreme photo-doping limit (n d = 0.5), the effective model (Ĥ dh;ex þĤ V ) becomes equivalent to the XXZ model 42 . Namely, we havê where J XY = − J ex , J Z = −J ex + 4V andη represents the ηoperators (see Methods). The XXZ model shows XY order (quasi-long-range order) for |J XY | > J Z , while it shows Ising order (long-range order) for |J XY | < J Z . In our language, the former corresponds to η-SC and the latter to CDW. Thus, it is natural that the phase transition between η-SC and CDW occurs at for strong photo-doping. Interestingly, the phase boundary remains located near this value over a wide photo-doping range, see Fig. 1.
As Fig. 2b-d shows, more than one order can be quasi-long ranged for a given set of parameters. When V is small, the SC correlations are dominant. The spin correlations are also quasilong-ranged, while the charge correlations show no sign of CDW (alternation of signs). When V is increased, the exponent of the spin correlations remains almost unchanged, while the decay of the SC correlations becomes faster, and the CDW correlation starts to develop around V ≳ 0.1 ¼ J ex 4 .
is in a coexistence regime where CDW, SDW, and η-SC orders are simultaneously quasi-long ranged. While the precise boundaries of the coexistence regime are difficult to determine (see Supplementary Note 3), by V = 0.4 (= J ex ), the CDW correlations become dominant and the SC correlations decay exponentially.
Origin and properties of the photo-doped phases. The photodoped states exhibit unique properties. Firstly, the η-SC is absent in equilibrium, since in chemically-doped systems either The critical exponent is extracted by fitting the correlation functions (χ(r)) with C 1 =r 2 þ C 2 cosðqrÞ=r a , where q = 2n d π, q = (1 − 2n d )π, and q = π for charge, spin, and SC correlations, respectively. We use the spatial distance r ∈ [6,30] for the fitting range. The phase boundary is only schematic and a guide to the eye except for n d = 0.5 and V = J ex /2, where η-SC and CDW are degenerated (η-SC/CDW).
doublons or holons are introduced and hence χ sc (r) vanishes. On the other hand, one expects that even in chemically-doped states, CDWs can develop due to the instability of the Fermi surface. Figure 3a, b however show that the CDW correlations are much stronger in photo-doped than in chemically-doped states. Furthermore, the CDW correlations show incommensurate oscillations with q = 2n d π, see Fig. 2a. This indicates that holons and doublons do not bind in pairs. Instead, the holons (doublons) are located in the middle of neighboring doublons (holons), even though the doublon-holon interaction V dh J ex 4 À V is attractive, seeĤ dh;ex þĤ V . The absence of binding is also directly confirmed by the evaluation of the doublon-holon correlations (Supplementary Note 2). This situation is in stark contrast with semiconductors, where the attractive interaction between photodoped holes and electrons leads to condensation of electron-hole pairs (excitons) at low temperatures [20][21][22]47 .
Let us discuss in more detail the physical origin of the CDW phase. In the extended Hubbard model the interaction between doublons (and between holons) is V dd À J ex 4 þ V, and thus V dd = −V dh . To investigate how the interactions among doublons and holons affect the CDW formation, we artificially add an interaction between neighboring doublons and holons, H V dh ΔV dh ∑ i ½n i;dniþ1;h þn i;hniþ1;d , so that the doublon-holon interaction becomesṼ dh V dh þ ΔV dh . We choose parameters such that bothṼ dh and V dd are repulsive. Figure 3c shows that the relative magnitude of the doublon-doublon (holon-holon) and doublon-holon interaction controls the physics and that an attractive doublon-holon interaction is not essential for the CDW. Namely, oscillations in χ c appear if V dd ≳Ṽ dh . This indicates that CDW correlations develop between the doublons and holons as if no singlons existed between them. The situation is analogous to the spin correlations in the one-dimensional t-J model, which can be explained by the squeezed Heisenberg chain without doublons and holons 59,60 . Underlying this phenomenon in the one-dimensional t-J model is the conservation of the spin configuration in the J → 0 limit 59 . Since a singlon always encounters the same neighbors, the system favors the spin configurations described by the Heisenberg hamiltonian. The same situation is realized in the photo-doped case. In the limit of J ex → 0, the configuration of doublons and holons is also conserved due to their peculiar kinematics, seê H kin;holon þĤ kin;doub . (Note that in normal semiconductors, holes and electrons can switch positions even in the one-dimensional case.) Thus, the configurations of doublons and holons are determined by the interaction termĤ dh;ex þĤ V , as in the case of spin configurations in the t-J model.
To confirm the above scenario, we evaluate the correlations between the doublons and holons in terms of a reduced distance which ignores singlons. The corresponding correlation function is defined asχ η z ðrÞ ¼ h∑ l ≥ 0 Q lη z rþlη z 0 Q l i 41 , where Q l is the Fig. 2 Correlation functions of the photo-doped states. Correlation functions for charge (χ c ), spin (χ s ), and superconductivity (χ sc ) are evaluated with the infinite time-evolving block decimation for the photo-doped states described by the effective HamiltonianĤ eff2 . a Normal scale plot for the nonlocal interaction V = 0.2 and the corresponding fitting with C 1 =r 2 þ C 2 cosðqrÞ=r a , where q = 2n d π, q = (1 − 2n d )π, and q = π for charge, spin and SC correlations, respectively. Here, r is the spatial distance, n d is the doublon density, and C 1 , C 2 , and the critical exponent a are fitting parameters. b-d Log-scale plots of the absolute value of the correlation functions for specified values of V. Empty (filled) markers correspond to χ < 0 (χ > 0). The dashed (dot-dashed) lines show C 2 r −a (C 1 r −2 ) extracted by fitting. For all panels, we use n d = 0.23 and r ∈ [6,30] for the fitting range. projection to states with l singlons between the 0th site and the (r + l)th site andη z i ¼ 1 2 ðn i;d Àn i;h Þ [ Fig. 3d]. Staggered correlations appear for V dd ≳Ṽ dh , which supports the above argument. Note that, even when doublons and holons show Ising-type order in the squeezed space, the correlations can still exhibit a powerlaw 61 . We thus conclude that the photo-induced CDW originates from the less repulsive doublon-holon interaction (compared to interactions between the same species), the peculiar kinematics of carriers and the one-dimensional configuration. Furthermore, the development of correlations between the doublons and holons in the squeezed system without singlons should also apply to systems with V ≲ J ex 2 . In these cases, the X and Y components of H dh;ex þĤ V (we regardĤ dh;ex þĤ V as an XXZ model, as in the extreme photo-doing limit) is dominant and the η-pairing phase emerges. This naturally explains the observation that the boundary between the η-pairing phase and the CDW phase is close to V ' J ex 2 independent of the photo-doping level.
Single-particle spectra. We now focus on the single-particle spectra, to clarify characteristic features of the different phases. Figure 4 shows the momentum-integrated spectrum A loc (ω) and the momentum-resolved spectrum A k (ω) for the η-SC phase [ Fig. 4a, b] and the CDW phase [ Fig. 4c, d] 46 . For the CDW phase, we use V = 1 to enhance the characteristic features. Unlike in equilibrium, but similar to photo-doped semiconductors, the photo-doped system exhibits two "Fermi levels" separating occupied (electron removal spectrum) from unoccupied (electron addition spectrum) states [ Fig. 4a, c]. The occupied states in the upper Hubbard band (UHB) region correspond to the removal of a doublon already existing in the photo-doped state, while those in the lower Hubbard band (LHB) region correspond to adding a holon (see Methods for precise definitions). In the η-SC phase, within our numerical accuracy, no gap signature appears in A k (ω) around the new Fermi levels [ Fig. 4b], which is in stark contrast to a normal superconductor with a gap around the Fermi level. The absence of a gap is also found for η-pairing states in higher dimensions 62 , and this suggests that the η-SC state is a kind of gapless superconductivity. On the other hand, in the CDW phase, gaps appear at the new Fermi levels [ Fig. 4d], as in the excitonic phase in photodoped semiconductors 22 . Finally, we observe that in-gap states between the UHB and LHB develop with photo-doping, which are more prominent for larger V [Fig. 4c]. These states may enable recombination processes suggesting that our assumption of approximately conserved doublon and holon numbers may become less valid as the excitation density increases. However, one needs to keep in mind the following points: (i) For large enough U, the Mott gap remains clear and the doublon and holon numbers are approximately conserved. Since the CDW is driven by V, the value of U does not affect its existence and spectral features. (ii) Even when in-gap states develop, an effective equilibrium description is meaningful. The recombination rate for a given state can be estimated by Fermi's golden rule, and if this rate is small compared to the intraband relaxation, the transient state can be described by (time-dependent) effective temperatures and chemical potentials 63 . Hence, our results show that the effective equilibrium description can be useful to study the closure or shrinking of a Mott gap via photo-doping, as a result of screened interactions and photo-induced spectral features 64 .

Discussion
We introduced a GGE-type effective equilibrium description for photo-doped strongly correlated systems. This provides a theoretical framework for systematic studies of nonthermal phases. Using this effective equilibrium description, we revealed emerging phases in the photo-doped one-dimensional extended Hubbard model. The η-pairing phase is stabilized in the small V regime even when the SU c (2) symmetry that protects η-pairing in the pure Hubbard model is absent, and it is characterized by gapless spectra. The CDW phase emerges in the larger V regime, and it is characterized by gapped spectra. These states are unique to a photo-doped strongly correlated system, where the peculiar kinematics of doublons and holons stabilizes them in a wide doping range. The similarity between the GGE-type description for strongly correlated systems and the pseudoequilibrium description for the photo-doped semiconductors allowed us to clarify some fundamental differences between these two systems. In particular, our results demonstrate that photo-doped strongly correlated systems and semiconductors exhibit qualitatively different phases due to the different nature of the injected carriers. Target systems to look for the characteristic Mott features include candidate materials of one-dimensional Mott insulators ranging from organic crystals, e.g., ET-F 2 TCNQ, to cuprates, e.g., Sr 2 CuO 3 , as well as cold-atom systems.
We also note that further insights into our results may be obtained by using equilibrium concepts such as Luttinger liquid theory and the exact wave function for J ex → 0 59 . In the nonequilibrium state, we expect at most three degrees of freedom: spin, pseudo-spin (consisting of doublon and holon), and charge (position of singlons). Thus, unlike in the equilibrium case, the maximum value of the conformal charge c would be 3, which may be realized in the η-pairing phase. A systematic analysis in this direction is under consideration.
The GGE-type description of photo-doped Mott states can be applied to various models and implemented with different equilibrium techniques such as slave boson approaches 65,66 and variational methods 51,67,68 . It can provide useful insights into experimental findings, e.g., photo-induced SC-like state, as well as theoretical results. For example, a recently found metastable orbital order in a photo-doped multi-orbital system can be reasonably explained by an effective equilibrium picture 69 . Systematic explorations of nonequilibrium phases in higher dimensions and at nonzero effective temperatures with the GGE-type description should be undertaken. For example, the GGE-type description allows us to investigate important aspects of nonequilibrium states, such as the screening of the interactions 64 or the stability of photoinduced phases against light irradiation 70 . These questions are interesting topics for future investigation.

Methods
Infinite time-evolving block decimation. The infinite time-evolving block decimation (iTEBD) method expresses the wave function of the system as a matrix product state (MPS), assuming translational invariance 51 . iTEBD directly treats the thermodynamic limit and we use cut-off dimensions D = 1000-3000 for the MPS to get converged results. We use the conservation laws for the numbers of spin-up and spin-down electrons at half-filling to improve the efficiency of the calculations. Away from half-filling, the trick with the conservation laws cannot be used, which makes the simulations less efficient. Thus, we use the exact diagonalization method in Fig. 3.
We also note that, to fit correlation functions obtained from iTEBD, we use functions with a power law. To be strict, because of the SU(2) spin symmetry, there may be a logarithmic correction on top of the power-law decay for the spin correlation. Still, in our results, fits with and without this correction work equally well. The latter fit yields larger a, hence the phase diagram and the qualitative arguments given in the main text remain unaffected.
Single-particle spectrum. The single-particle spectrum is defined as follows. The local spectrum is A loc ðωÞ À 1 π ImG R i ðωÞ and the momentum-resolved spectrum is A k ðωÞ À 1 π ImG R k ðωÞ. Here, G R i ðωÞ (G R k ðωÞ) is the Fourier transform of the retarded Green's function G R i ðtÞ ¼ Àih½ĉ iσ ðtÞ;ĉ y iσ ð0Þ þ i (G R k ðtÞ ¼ Àih½ĉ kσ ðtÞ;ĉ y kσ ð0Þ þ i). Note thatĉ iσ ðtÞ is the Heisenberg representation ofĉ in terms ofĤ eff . The occupied spectra correspond to A < loc ðωÞ 1 2π ImG < i ðωÞ and A < k ðωÞ 1 2π ImG < k ðωÞ, where G < denotes the lesser part of the Green's functions. To evaluate these quantities using the effective equilibrium description and iTEBD, we employ the method proposed by some of the authors 46 . Namely, we evaluate G R (t) using an auxiliary band and perform the Fourier transformation with a Gaussian window, F Gauss ðtÞ ¼ exp À t 2 2σ 2 À Á , and we use σ = 5.0. Thus, the broadening of the resultant spectrum is inevitable and a gap much smaller than the broadening cannot be captured. Still, we checked that no gap signature appears in the η pairing state for an increased value of J ex , where we would expect an increase of the gap (if any).
Here, we introduce the η-operators asη þ i ¼ θ iĉ y i#ĉ y i" ,η À i ¼ θ iĉi"ĉi# and η z i ¼ 1 2 ðn i À 1Þ with θ i = (−) i . The shift of the local interaction is described bŷ H Here, the superscript "(2)" indicates that the term is order of O Here, 〈k, i, j〉 means that both of (k, i) and (i, j) are pairs of neighboring sites. The sum is over all possible such combinations (without double counting), where we regard 〈k, i, j〉 = 〈j, i, k〉.
In the evaluation of the physical quantities, we use the operators of the effective model. To be strict, if physical quantities for the original Hamiltonian are to be computed, one also needs to take account of corrections from the SW transformation to the operators. However, these corrections are not necessary to see the leading behavior. The same strategy is often used in the evaluation of physical quantities for the Heisenberg model or the t-J model.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Code availability