Enhancement of Local Pairing Correlations in Periodically Driven Mott Insulators

We investigate a model for a Mott insulator in presence of a time-periodic modulated interaction and a coupling to a thermal reservoir. The combination of drive and dissipation leads to non-equilibrium steady states with a large number of doublon excitations, well above the maximum thermal-equilibrium value. We interpret this effect as an enhancement of local pairing correlations, providing analytical arguments based on a Floquet Hamiltonian. Remarkably, this Hamiltonian shows a tendency to develop long-range staggered superconducting correlations. This suggests the possibility of realizing the elusive eta-pairing phase in driven-dissipative Mott Insulators.

We investigate a model for a Mott insulator in presence of a time-periodic modulated interaction and a coupling to a thermal reservoir. The combination of drive and dissipation leads to nonequilibrium steady states with a large number of doublon excitations, well above the maximum thermal-equilibrium value. We interpret this effect as an enhancement of local pairing correlations, providing analytical arguments based on a Floquet Hamiltonian. Remarkably, this Hamiltonian shows a tendency to develop long-range staggered superconducting correlations. This suggests the possibility of realizing the elusive eta-pairing phase in driven-dissipative Mott Insulators.
The Floquet engineering of complex quantum systems is a very active line of research in today's condensed matter physics [1]. It consists in the design of periodic perturbations to achieve non-equilibrium driven states remarkably different from their undriven counterparts. Examples are the dynamical control of band topology [2,3] and of magnetic interactions [4] in ultracold atoms in optical lattices, and of effective Hamiltonian parameters in solids under intense laser-pulse excitation [5][6][7].
A useful description of a periodically driven quantum system is in terms of effective static Hamiltonians derived with large-frequency expansions [8][9][10][11]. In general, however, the drive affects also the distribution function of the system, eventually leading to thermalization to a trivial infinite-temperature state [12,13]. Nevertheless, when heating can be avoided for finite but long times, interesting prethermal Floquet states can be observed. This is the case, for example, for very large drive frequency [14][15][16][17] or systems close to integrability [18][19][20][21][22][23]. In particular, Ref. [24] showed that strong electronic correlations lead to finite-frequency prethermal states with remarkable properties as a function of drive frequency.
A natural question concerning the Floquet prethermal state is whether the coupling to external reservoirs would cancel out its interesting features, or rather preserve them and possibly make them more accessible. Particularly interesting is the possibility to control the distribution function of the system by means of a dissipation mechanism of the energy injected by the drive [25][26][27][28][29].
To investigate this point, in this work we consider the Fermi-Hubbard model with a periodically driven interaction and coupled to a thermal reservoir. Starting from the Mott-insulating phase, our numerical calculations show that the combination of drive and dissipation leads to steady states that are not accessible in the corresponding isolated model. In particular, we reveal a regime with a remarkably large number of high-energy doublon excitations, well above the maximum equilib-rium value for the half-filled repulsive Hubbard model.
We interpret this steady-state large double occupancy as an enhancement of local pairing correlations, and we describe the effect as a thermalization to a lowest-order Floquet Hamiltonian. Remarkably, we find that higherorder terms promote finite-momentum doublon superfluidity, namely staggered long-range pairing correlations among fermions (η-pairing), which spontaneously break the hidden SU C (2) charge symmetry of the half-filled Hubbard model [30][31][32]. This suggests a nonequilibrium protocol for Floquet engineering exotic superconducting states in driven-dissipative Mott insulators, as also very recently investigated in similar contexts [33][34][35][36].
Our results are relevant for current experiments on laser-pumped organic Mott insulators [5,6] and ultracold Fermi gases in driven optical lattices [37,38]. We discuss the latter in particular, suggesting to explore a possibly overlooked regime in future experiments.
Model -The Hamiltonian of the driven-dissipative Fermi-Hubbard model reads H = H Hub + H diss , where: Here the c's operators describe fermions hopping with amplitude V ij and subject to a driven local interaction U (t ≥ 0) = U 0 + δU sin Ωt. The bare density of states is semicircular with bandwidth 4V and we measure energy, frequency and inverse of time ( = 1) in units of V [39]. The thermal bath is implemented by independent sets of bosonic modes b's which couple to density at each lattice site, with spectral function J(ω) = α g 2 α δ(ω − ω α ) ∝ ω 2 e − ω ωc (ω c = 1) and coupling λ [39]. Importantly, the bath allows energy dissipation but commutes with density n i = n i↑ + n i↓ and preserves particle-hole symmetry. The system remains half filled at all times ( n iσ = 0.5). Starting from a thermal-equilibrium state, we calculate the time evolution by means of nonequilibrium dynamical mean-field theory (DMFT) [40,41] with the noncrossing approximation as impurity solver [42,43], including the effect of dissipation at weak coupling in λ [44] (see Supplemental Material [45] Sec. I for implementation details, and Sec. II for one-crossing-approximation benchmarks). We calculate double occupancy D(t) = n i↑ (t)n i↓ (t) , kinetic energy K(t) = V ij c † iσ (t)c jσ (t) , and local Green's function G σ (t, t ) = −i Tc iσ (t)c † iσ (t ) . For definiteness, we choose U 0 = 8 for the initial Mottinsulating state in equilibrium at T = 1 [46] and δU = 2 for the drive amplitude (see Supplemental Material [45] Sec. VI for a discussion on δU ). The bath temperature is T bath = 1 unless specified differently. In absence of dissipation, Floquet prethermalization is observed at all frequencies except for the resonance Ω * = 8.12 U 0 [24]. We restrict ourselves to paramagnetic states, leaving the interplay of drive, dissipation and magnetism to future studies.
Time evolution -In the driven-dissipative model, as well as in the isolated case, double occupancy and kinetic energy display a separation of time scales between fast oscillations synchronized with the drive and a slowly varying average value. However, after a common transient, the thermal reservoir starts to be effective and changes substantially the long time behavior of both observables.
For weak bath coupling and drive above resonance, the double occupancy grows substantially larger than in the isolated model, going to a stationary average above 0.25 ( Fig. 1a, λ = 0.2). Such a large value would be possible, at equilibrium, only if the interaction were attractive. This striking effect highlights the peculiarity of this nonequilibrium steady state, as we discuss thoroughly below.
Upon increasing the bath coupling (Fig. 1a, λ = 1.0), the double occupancy decreases and eventually remains below the limit of 0.25 at all times. Moreover, we notice that the bath is effective only after a transient time ∼ 1/λ 2 , which makes the regime of very weak coupling not accessible by the numerical simulation (see also Ref. [47]).
At the same time, the kinetic energy is also largely affected by dissipation (Fig. 1b). Here the effect is more intuitive: in the isolated model the drive leads to a prethermal state with positive kinetic energy, indicative of a population inversion [22,24]. On the other hand, the thermal reservoir dissipates the excess kinetic energy, which remains negative as at equilibrium, and inhibits the population inversion, as we also explicitly show below.
Long-time average -To study the role of drive frequency and bath coupling in a more systematic way, we consider the long-time average of double occupancy and kinetic energy. For weak bath coupling, the dissipative model has double occupancy larger than the isolated one at all frequencies (Fig. 1c, λ = 0.2). However, a remarkable change happens crossing the resonance Ω * U 0 . Below resonance, the dissipation has only a quantitative, rather weak effect. In contrast, above resonance, we systematically observe a large increase of double occupancy across the limit of 0.25, as discussed previously for a selected frequency. Lower values are then recovered upon increasing the frequency further, as the system eventually becomes transparent to the drive.
Independent of the bath coupling, the kinetic energy of the dissipative model is rather featureless and negative for all frequencies (Fig. 1d, λ = 0.2, 1.0). Thus, the thermal reservoir cancels the region of positive kinetic energy characteristic of the isolated case ( Fig. 1d, λ = 0.0).
The difference between below and above resonance appears also in the dependence on bath coupling (Fig. 1e) and bath temperature (Fig. 1f). Below resonance (Ω = 7) the double occupancy is almost independent of bath coupling and decreases on lowering the bath temperature. Quite differently, above resonance (Ω = 9) it increases on lowering the bath temperature at weak bath coupling.
Notice that the observed behavior does not depend on the details of the bath spectral function as long as we deal with bosonic modes [39]. In contrast, a fermionic reservoir does not lead to the same steady-state large double occupancy (see Supplemental Material [45] Sec. III) because it can change the local density even at zero hopping, spoiling the quasi-conservation of doublons which is the basis of Floquet prethermalization in this system [24].
Spectral function -To gain insight into the nature of the steady state, we calculate the spectral functionĀ(ω) and occupation functionN (ω) as the average Wigner transforms of the retarded and lesser components of the local Green's function [24]. While the spectral function is the same in the isolated and dissipative models, the occupation function, and thus the distribution function F (ω) =N (ω)/Ā(ω), changes drastically for drive frequency above resonance and weak bath coupling.
In the isolated model,N (ω) is shifted towards high energy with respect toĀ(ω) (Fig. 2a). There is therefore a population inversion within the Hubbard bands. Indeed, the local behavior ofF (ω) for ω ±U 0 /2 has the shape of a Fermi function with negative temperature (Fig. 2c).
The thermal reservoir completely changes the situation. First, as the dissipation enhances the energy redistribution within the Hubbard bands,N (ω) is pushed back to lower energy (Fig. 2b), cancelling the population inversion. As a consequence,F (ω) assumes the shape of a Fermi function with positive temperature for ω ±U 0 /2 (Fig. 2d). Then, the overall weight ofN (ω) in the upper band grows and becomes even larger than in the lower band, meaning the creation of a large number of highenergy doublon excitations. These two effects are qualitatively related to the ones discussed above: change of sign of kinetic energy and growth of double occupancy.
Discussion -The above numerical results demonstrate that, in the strongly repulsive Fermi-Hubbard model, the combination of a time-periodic interaction and a dissipative bath leads to steady states with a remarkably large number of doublon excitations. Interestingly, this large double occupancy immediately translates into enhanced local pairing correlations D = c † i↑ c † i↓ c j↓ c j↑ | i=j , although a full calculation of the lattice susceptibility within DMFT is beyond the scope of this paper.
In order to unveil the origin of this effect, let us first consider the isolated model and its Floquet Hamiltonian. To this end, we consider a frequency close to resonance Ω Ω * U 0 and perform a rotating-frame transformation on the Hamiltonian (1), followed by a high-frequency expansion [9,48] (see Supplemental Material [45] Sec. V). At lowest order we find: are those hopping terms in Eq. (1) that do not alter the number of doubly occupied sites (n iσ = 1 − n iσ ,↑ =↓ and↓ =↑). Eq. (3) can be interpreted as a Hamiltonian of doublons and holons, where the first term contains hopping and the second term acts as a chemical potential. The numerical results around Ω * U 0 in the isolated model are qualitatively captured in terms of thermalization to this effective Hamiltonian. Indeed, we can extract the effective temperature T eff (Ω) (U 0 −Ω) −1 (see Supplemental Material [45] Fig. S3) and, since |T eff | V for Ω Ω * , we can disregard the kinetic term in Eq. We now turn to the dissipative model. Here two observations are crucial. First, the enhancement of double occupancy is most pronounced for weak bath coupling (see Fig. 1e). Second, the dissipation leaves largely unchanged the spectral function of the system, while it profoundly changes its occupation (see Fig. 2). On this basis, we argue that at weak coupling the bath does not change the Floquet Hamiltonian, but only affects the effective temperature T eff T bath = 1 (see Supplemental Material [45] Fig. S3). This leads to D 0.25+(Ω−U 0 ) which qualitatively reproduces the behavior for λ = 0.2 around Ω Ω * (Fig. 1c; dashed line) where again the mismatch is due to the finite hopping in Eq. (3).
The outcome of this analysis is that, at least for moderately high temperatures T bath V , the driven-dissipative protocol of Eqs. (1) and (2) leads at long times to thermal states of the doublon Hamiltonian (3). Singly occupied sites, which are relevant in the transient dynamics during doublon-holon proliferation, are not relevant for the steady-state physics, and can be considered as a reservoir of energy and particles to the doublon-holon system.
A natural question now is whether the enhanced local pairing correlations can propagate through the lattice giving a superfluid state of doublons. To answer this, we consider the next order in the Floquet Hamiltonian [45]: (4) Here J n (x) is the n-th order Bessel function of the first kind, K + = (V ij /V )c † iσ c jσ n iσnjσ = (K − ) † and one has to note that, in the case of weak drive amplitude considered here, J n (δU/Ω) ∼ (δU/Ω) n and therefore all terms in Eq. (4) indeed vanish as the inverse drive frequency.
The first two terms in parentheses in Eq. (4) create or annihilate doublon excitations, controlling the transient dynamics. However, these processes are largely inhibited in the steady state. Indeed, these terms depend strongly on the drive amplitude δU , which controls the transient time-scale but does not influence the long-time steadystate, as found in both the isolated [24] and the dissipative models (see Supplementary Material [45] Sec. V).
The last term in Eq. (4) is similar to the Schrieffer-Wolff result [48,49], which is retrieved for δU → 0 at fixed Ω = U 0 . It contains several contributions such as density and exchange interactions, and correlated threesite hopping processes. At equilibrium and half filling, in the relevant limit of zero double occupancy, this gives the usual anti-ferromagnetic Heisenberg model. In farfrom-equilibrium situations, if a sizeable population of doublons is achieved as in the present case, the same term leads to a completely different physics (see also Ref. [50]).
To discuss Hamiltonian (4) on states with very large double occupancy, it is instructive to neglect processes involving singly occupied sites and rewrite it as [45]: Here J eff = 2V 2 /Ω(J 0 (δU/Ω)) 2 , the first term is a doublon hopping, and the second term a first-neighbor doublon interaction. It is now convenient to consider a transformation on spin-down operators c i↓ →c i↓ = (−1) i c † i↓ which recasts Eq. (5) asH eff Hub = −J eff η i ·η j namely as an isotropic ferromagnetic Heisenberg model for the socalled η-spins η i = 1 2 αβc † iα σ αβciβ [45]. The invariance under η-spin rotation is associated to the charge SU C (2) symmetry of Hamiltonian (1) which can be used to build eigenstates of the Hubbard model with staggered longrange superconducting correlations (η-pairing) [30,31]. These same correlations are encoded in Eq. (5). Indeed, the η-spin ferromagnetic Heisenberg model has magnetization η z = D − 0.5 and below a critical temperature ∼ J eff it has finite order parameter in the xy-plane, which corresponds to staggered long-range pairing correlations . We stress again that here, for simplicity, we do not consider the interplay between doublons and singly occupied sites, which would result in additional terms in Eq. (5). Moreover, we notice that the SU C (2) symmetry implies a degeneracy between the xy-plane and the z-axis of the η-spin, which translates into a competition between superfluidity and charge-density wave. We leave the investigation of these issues for future work.
The model system investigated here can be realized in current experimental platforms. Particularly promising are Mott-insulating organic molecular crystals, where laser excitations can induce an effective time-periodic modulation of interaction [5,6]. More direct control is achieved with ultracold atoms in optical lattices. Recent experiments [37,38] have studied the Floquet prethermal state and remarkably found large double occupancy for drive above resonance. We suggest that also in these cases a key role is played by dissipation, which is unavoidable even in cold atoms. Finally, we notice that these experiments have focused on the regime where the effective Hamiltonian reduces to a renormalized Hubbard model. In contrast, here we have studied the case of a doublon-only Floquet Hamiltonian. Therefore, we suggest future experiments to investigate this latter regime to detect the presence of staggered pairing correlations.
Conclusions -In this work, we have studied the combination of a periodically driven interaction and a dissipative bath in the strongly repulsive Fermi-Hubbard model. For weak bath coupling and frequency in a range above the resonance of the isolated model [24], we find a large increase of double occupancy, well above the maximum equilibrium value, which we interpret as an enhancement of local pairing correlations, and understand in terms of thermalization to the lowest-order Floquet Hamiltonian.
Remarkably, the next-order Floquet Hamiltonian contains terms which promote staggered pairing correlations. Therefore, provided a nonequilibrium protocol to reach low enough effective temperatures (see e.g. Ref. [51]) and eventually further increase the doublon density, the steady state of the driven-dissipative Fermi-Hubbard model would contain off-diagonal staggered long-range order (η-pairing), hence a superfluid phase of doublon excitations, similarly to what very recently found in similar models [35,36]. Here we describe our implementation of a bosonic bath in the non-crossing-approximation (NCA) impurity solver for non-equilibrium dynamical mean-field theory (DMFT). Starting from Hamiltonians (1) and (2) in main text, we integrate out the bosons and obtain the action: (S1) Here the integrals are along the three-branch Keldysh contour and the bath enters via the hybridization where θ(t, t ) is the Heaviside theta function on contour and n b (ω) is the Bose distribution.
In DMFT the lattice action (S1) is mapped onto the action of a quantum impurity coupled to a self-consistent fermionic bath: Here we have used the relation ∆(t, t ) = V 2 G(t, t ) for the hybridization ∆(t, t ) of the self-consistent fermionic bath, which is valid on Bethe lattice. To derive the NCA equations, we expand the partition function Tr(exp[−iS imp ]) into a power series in V and λ and truncate the expansion at the first self-consistent order [42][43][44]. This series is expressed in terms of the propagator of the states of the impurity R, and of its self-energy S, which satisfy an integro-differential equation similar to the usual Dyson equation (see Ref. [24] for our implementation). Then, in the present case a convenient basis choice is {|0 , |↑ , |↓ , |↑↓ } which makes the propagator R and the self-energy S diagonal. Finally, we exploit the spin SU(2) symmetry and the particle-hole SU C (2) symmetry of Eq. (S2) which further reduce the number of propagators to two: one for the empty state |0 and one for the singly-occupied state |↑ , with the following self-energies:

II. OCA BENCHMARK
In this work we investigate the Mott-insulating phase of the Fermi-Hubbard model with average interaction parameter U 0 = 8V much larger than the bare band-width W = 4V . Moreover, we consider initial thermal density matrix at a rather high temperature T = 1. In these conditions, the NCA approximation is expected to perform well both at equilibrium and out of equilibrium [42].
To confirm this, we have performed some calculations using the next-order one-crossing approximation (OCA). This takes into consideration terms of order O(V 4 ) in the hybridization expansion [42,43] (see Ref. [24] for our implementation). As expected, the time evolution of double occupancy and kinetic energy (Fig. S1) is essentially identical to the one shown in the main text ( Fig. 1a-b).

III. FERMIONIC BATH
It is interesting to consider an external reservoir of fermionic modes, instead of the bosonic bath considered in the main text. To do so, we substitute Eq. (2) of the main text with the following coupling to a fermionic bath: Here the operators f 's represent sets of independent fermionic harmonic oscillators. In this case, the system can dissipate both energy and particles. Moreover, this type of coupling can be treated exactly: integrating out the fermionic bath, we introduce an additional hybridization ∆ f (t, t ) which can simply be added to the DMFT self-consistent hybridization ∆(t, t ). As shown in Fig. S2, the coupling with a fermionic bath does not lead to the same interesting features discussed in the main text for the bosonic bath (cf. Fig. 1a-c). This is due to the fact that the bosonic bath, as opposed to the fermionic bath, commutes with the local density and conserves double occupancy, and therefore it preserves the mechanism for Floquet prethermalization.

IV. EFFECTIVE TEMPERATURE
In Fig. 2 of main text we show the steady-state average distribution functionF (ω) for Ω = 9 > Ω * , along with a Fermi-function fit [1+exp(ω/T eff )] −1 around the center of the upper Hubbard band ω U 0 /2. From this fit we can extract the effective temperature T eff . Here, in Fig. S2 we plot this effective temperature for the isolated and the dissipative models, as a function of the drive frequency.

V. ROTATING-FRAME TRANSFORMATION AND HIGH-FREQUENCY EXPANSION
To carry out the large-frequency expansion [9,48], we first need to transform Hamiltonian (1) of main text to a rotating frame with respect to the interaction: F (t) = Ωt − (δU/Ω) cos Ωt.
It is useful to introduce the operators K 0 and K ± as in the main text (n iσ = 1 − n iσ ,↑ =↓ and↓ =↑): In the former, the effective temperature diverges as (Ω * − Ω) −1 (dashed line). The divergence at Ω = Ω * signals the resonant thermalization, while the negative effective temperature for Ω > Ω * signals the population inversion [24]. In contrast, the dissipative model has temperature T bath (solid line) fixed by the thermal reservoir.
The Fourier components H m are given in Eqs. (S15) and depend on frequency through the Bessel function. When the large-frequency limit is taken at fixed δU/Ω, this dependence does not show up in the expansion. In contrast, here we keep δU constant and we have to consider the asymptotic behavior J n (δU/Ω) ∼ (δU/Ω) n . Then, there are terms in the average Hamiltonian (S15a) which vanish as Ω −1 and do not enter the lowest order of the expansion (cf. Eq. (3) of main text): For the same reason, at first order only enter those terms which actually vanish as Ω −1 (cf. Eq. (4) of main text): (S18)

Eq. (5) of main text
Here we calculate the commutator in Eq. (S18): (S19) The commutators in Eq. (S19) give three-and two-site terms. If we retain only the two-site terms, then the first sum in Eq. (S19) reads: (S20) Now, with the identities c iσ n iσ = c iσ and n iσ c iσ = 0, together with their Hermitian conjugates, it is easy to see that the terms in the second sum in Eq. (S19) are non-vanishing only if {i = l, j = k, σ =σ}, giving:

VI. ROLE OF DRIVE AMPLITUDE
In Fig. S4 we compare the time evolutions for various drive amplitudes δU . These show interesting features in both the isolated model (see also Ref. [24]) and the dissipative model. First, it is evident how the drive amplitude controls the time scale of the transient dynamics, with larger amplitudes inducing shorter transients. Second, the steady steady, in sharp contrast, appears to be largely independent from the drive amplitude.
These observations are consistent with and provide further validation to our analysis based on the Floquet theory (Eqs. (3) to (5) of main text). Indeed, in Eq. (4) a factor J 1 (δU/Ω) ∼ δU/Ω multiplies the terms creating doublon excitations, which are responsible for the increase of double occupancy in the transient dynamics. Thus, a large drive amplitude enhances these terms and makes the transient shorter. Moreover, since the long-time values do not depend on δU , it is confirmed that these terms are largely suppressed in the steady state, which validates the introduction of the doublon-only Hamiltonian (5).

VII. HYPERCUBIC LATTICE
In the main text we consider a system on Bethe lattice, whose free-electron density of states is semicircular. Here we consider a hypercubic lattice in infinite dimensions, whose free-electron density of states is Gaussian. In Fig. S5 we show that this difference does not qualita-tively affect the physics discussed in the main text. For drive frequency above resonance Ω = 9 > Ω * U 0 we observe qualitatively the same increase in the long-time double occupancy which goes above 0.25 in the dissipative model (λ = 0.2) both on Bethe lattice (Fig. S5a) and on hypercubic lattice (Fig. S5b). The spectral and occupation function on hypercubic lattice also show the same changes between the isolated (Fig. S5c) and the dissipative (Fig. S5d) models. In the former, we observe population inversion within each of the Hubbard bands. In the latter, as a consequence of dissipation, the population inversion within the Hubbard bands does not take place. At the same time, the occupation of the upper Hubbard band grows larger than the one of the lower Hubbard band (cf. Fig. 2 of main text).

VIII. BATH SPECTRAL FUNCTION
In the main text we consider a bosonic thermal bath with spectral function J(ω) = ωc 2 ( ω ωc ) s e − ω ωc with ω c = 1 and s = 2. Here in Fig. S6 we provide additional numerical results with various parameters ω c and s.