Heating-Induced Long-Range $\eta$-Pairing in the Hubbard Model

We show how heating the spin degrees of freedom of the Hubbard model to infinite temperature can be used to melt the order within this sector and reach steady states, in any dimension, which have completely uniform long-range correlations between $\eta$ pairs. We induce this heating with either dissipation or periodic driving and evolve the system towards a non-equilibrium steady state. This steady state is identical in both cases and displays distance-invariant off-diagonal $\eta$ correlations. These correlations were first recognised in the superconducting eigenstates described in C. N. Yang's seminal paper (Physical Review Letters, 63, 2144 (1989)), which are a subset of our steady states. Finally, we show that our results are a consequence of the symmetry properties of the model and independent of the microscopic details of the heating mechanism.

Introduction -Driving and dissipation have recently emerged as transformative tools for dynamically evolving quantum systems into non-equilibrium phases with desirable properties [1][2][3][4][5]. In the context of stronglycorrelated many-body systems these tools can drastically alter their microscopic behaviour and manifest a variety of collective and cooperative phenomena at the macroscopic level.
Controlled dissipation, for example, has been shown to be a versatile resource for quantum information and simulation purposes [6,7]. Recent proposals show how, in ultracold atomic systems, Markovian baths which act quasi-locally can drive the system towards pure steady states with exotic properties such as superfluidity [8,9]. These schemes, however, rely on precise engineering of the Lindblad jump-operators in order to target specific states and avoid the system heating up to a generic thermal ensemble.
Meanwhile, significant interest has been generated by the superconducting-like states that have been induced in a variety of materials by transient excitation of the vibrational degrees of freedom using THz laser pulses [10,11]. These optically driven systems are often modelled via Floquet driving -where the Hamiltonian is subject to a time-dependent periodic field [12][13][14]. Through careful choice of the driving parameters the effective Hamiltonian can be modified, on comparatively short timescales, to one which favours superconductivity and so the system may transiently reach a superconducting prethermal state. However, Floquet heating [15][16][17] means that in most cases these systems continue to absorb energy from the driving field, heating them up and eventually destroying any superconducting order present.
In this Letter we show how, counter-intuitively, the interplay between symmetry and heating can actually lead to persistent, coherent, long-range correlations. We achieve this heating either with dephasing or with periodic driving and the symmetry means the microscopic details of these mechanisms are unimportant. For concreteness, we apply them to the spin degrees of freedom of bi-partite D-dimensional realisations of the Hubbard model, melting any order in this sector and reaching robust mixed states with completely uniform long-range correlations between η-pairs. This pairing is known to provide a possible mechanism for superconductivity [18][19][20]. In the case of dephasing we prove the results of our simulations by using symmetry arguments to blockdiagonalise the Liouvillian superoperator, giving an explicit parametric form for the steady states.
Model -The Hubbard model, in second quantised form, reads where c † σ,i and its adjoint are the usual creation and annihilation operators for a fermion of spin σ ∈ {↑, ↓} on site i. Additionally n σ,i is the number operator for a particle of spin σ on site i and τ , U and µ play the role of kinetic, interaction and chemical energy scales respectively. Here, and in the remainder of this work, we have set = 1.
The Hamiltonian in Eq. (1) has a rich structure due to the two SU(2) symmetries it possesses [21]. The first of these is the well known 'spin' symmetry which relates to spinful particles (singlons) σ ∈ {↑, ↓}. The second, which we refer to as 'η-symmetry' and central to this letter, relates to spinless quasi-particles (doublons and holons) σ ∈ {↑↓, vac}. This η-symmetry can be revealed by introducing the associated operators with η + i and η − i describing the creation and annihilation of a doublon on site i. The direction of its momentum is dictated by the parity of the site. The operators in Eq. arXiv:1902.05012v1 [quant-ph] 13 Feb 2019 (2) fulfill the relations [H, η ± ] = ±2µη ± and [H, η + η − ] = [H, η z ] = 0 and they commute with all the generators of the spin symmetry.
The presence of η-pairing superconductivity in the Hubbard model is a phenomenon first recognised by C. N. Yang in his seminal paper [22]. There it was proved that pure states of the form |ψ ∝ (η + ) N |vac are eigenstates of H and possess off-diagonal long-range order (ODLRO) in the form of doublon-doublon correlations where ρ = |ψ ψ|. This relation provides a possible definition of superconductivity as a finite value of this quantity can be shown to imply both the Meissner effect and flux quantisation [18,23,24]. These states, however, are excited states of H and the long-range order they possess is not usually seen in physical states (ground states, thermal states etc.) of the model due to destructive interference from the short-range coherences involving spinful particles. We will now show that by driving the Hubbard model in the spin basis these are destroyed and we are able to consistently engineer states with long-range uniform correlations in η + i η − j . As our primary mechanism for achieving this we consider the Hubbard model immersed in an environment which induces local dephasing in the spin basis. This model is motivated by the spin fluctuation theory of superconductivity, where electrons pair due to their scattering on spin fluctuations [25][26][27]. Our dephasing mimics these scattering events and recent experiments have demonstrated that spin-photon interactions can be induced by placing spin-electron ensembles within Magnonic cavities [28,29].
Despite the 'toy' nature of our model the Fermi-Hubbard Hamiltonian can be accurately realised by loading ultracold fermionic atoms into optical lattices [30][31][32]. These quantum simulators offer precise experimental control over the microscopic details of the system and dephasing may be induced by immersion of the lattice into an atomic bath in order to engineer scattering processes between fermions and the background atoms [33,34].
Results -We couple the Hubbard Hamiltonian in Eq. (1), H, to an environment which applies spin dephasing on each site of an M -site lattice. The ensuing dynamics is modelled, under the Markov approximation, via the Lindblad master equation with Lindblad operators s z k on each site. The time evolution is encoded in the Liouvillian superoperator L. Because the Lindblad operators are restricted to the spin SU(2) symmetry sector we obtain [L k , η ± ] = [L k , η z ] = 0 ∀k. For the Lindblad equation any operator A which satisfies is a null eigenvector of L and the adjoint map L † , which means the associated observable A is conserved. Thus, the map in Eq. (4) conserves η + η − , η z and S z = i s z i . Through a power series expansion, any appropriately trace normalised function f (η + η − , η z , S z ) of these operators is a steady state of L, i.e. if ρ ss = f (η + η − , η z , S z ) then Lρ ss = 0. In the Supplemental Material we show that steady states of this form have Tr(ρη + i η − j ) = const ∀i, j (see also [35]). We emphasize that this expected steady state is realisable for any bipartite D dimensional realisation of our lattice model.
Here, we focus on 1D lattice realisations due to their numerical tractability. We perform calculations where the system is initialised in the ground state of the Hubbard Hamiltonian for some τ and U . We then evolve the system by solving the time-dependent master equation in Eq. (4) and quenching U . For these simulations we perform a 'quantum trajectories' approach [36] to solving Eq. (4) combined with DMRG (Density Matrix Renormalisation Group) [37] and TEBD (Time Evolving Block Decimation) [38] algorithms for finding the ground state and evolving the wavefunction of the system respectively. These simulations were performed with the use of the Tensor Network Theory (TNT) library [39]. Further numerical details can be found in the Supplemental Material.
In Fig. 1. we plot, for various distances j, the averaged quantity η + i η − i+j over time for a quench from the ground state of the Hubbard Hamiltonian. We also include the matrix of correlations of η + i η − j for the density matrix following the quench. For comparison, we show the case where the system is not coupled to the environment, γ = 0. When the environment is present there is a relaxation of the system to a steady state which possesses the order expressed in Eq. (3); the value of η + i η − j is finite and constant for all i = j. We observe that this is facilitated by a decrease in the short-range η-pairing correlations in order to allow the long-range ones to increase: the non-equilibrium dynamics involve a 'spreading' of the ηcorrelations over all length scales of the system which is necessary due to the conservation of η + η − . Plots of the corresponding doublon momentum distribution are discussed in the Supplemental Material.
We also show in the Supplemental Material how this ODLRO is observable even under perturbative and unwanted dephasing in the model. We do this by introducing 'charge-dephasing' using jump operators that don't solely act in the spin basis of the model and cause the steady state of the Liouvillian to contain no coherences. Nonetheless, by setting the strength of the chargedephasing to values around 1% of the spin-dephasing, we find there still exists an intermediate time-scale (around ∆tτ = 30 hopping times) in which the results of Fig. 1. can be observed. For perturbations such as this, a quantum Zeno effect [40][41][42] is apparent as the duration of the window in which the uniform long-range correlations exist can be directly controlled by the ratio of the couplings for the spin dephasing and the charge dephasing. In particular, we show that the length of this window can be extended by increasing the coupling of the spin dephasing whilst the amplitude of the charge dephasing, and the other parameters in the model, remain constant.
To obtain a more intuitive idea of the nature of the steady state ρ ss = f (η + η − , η z , S z ) we consider a physically motivated parametrisation of the steady state function where we have exchanged η z and S z with N ↑ = i n ↑,i and N ↓ = i n ↓,i as they can be expressed as linear combinations of these operators. Calculations on small lattices show that the parametrisation in Eq. (6)  merical simulations which reach the long-time limit in Eq. (4). Equation (6) describes a generalized grandcanonical-like equilibrium state (GCE) with the Lagrange multipliers µ 1 , µ 2 , µ 3 associated to each of the conserved quantities. Notably, the Lagrange multiplier for the Hamiltonian is 0, i.e. at infinite temperature. As is typical of states in this form [43][44][45] the Lagrange multipliers are independent of the quench parameters and are governed solely by the values of the conserved quantities η + η − , N ↑ and N ↓ associated with the initial state.
To emphasize the properties of the stationary state in Eq. (6) we confine ourselves to symmetric half-filling N ↑ = N ↓ = M/2 and plot, in Fig. 2., the projection of the density matrix in the spin (s) and η (d) degrees of freedom P s,d ρP s,d . Explicitly, these projectors remove any basis vectors which contain, respectively, doublons (σ ∈ {↑↓, vac}) and singlons (σ ∈ {↑, ↓}) on any lattice sites. We plot these projections following a quench from a thermal state of the Hubbard model. Initially (Figs. 2a. and 2b.), the system contains coherences in both these sectors which decay with distance -additionally (not pictured) the density matrix contains coherences between doublons and singlons. Following spin dephasing in the long-time limit any coherences involving singlons are destroyed resulting in an infinite temperature ensemble in the spin symmetry sector (Fig. 2c.). Only coherences involving doublons and holons remain (Fig. 2d). These are completely distance invariant, as described by the parametrisation in Eq. (6).
The strong symmetries [46] η + η − , N ↑ and N ↓ can be used to block-diagonalise the Liouvillian and hence the degeneracy of ρ ss is determined by the distinct combinations of the eigenvalues of these operators. The steady states are the states with maximum η order in each block which we quantify by the parameter | η + i η − j |, i = j for the amplitude of the uniform off-diagonal correlations. Notably, Yang's states |ψ ∝ (η + ) N |vac form a subset of these steady states. They are the only pure steady states as they exist in the decoherence free subspace [2] of Eq. (4).
The infinite temperature ensemble pictured in Fig. 2c. shows that the dephasing has continuously pumped energy into the spin degrees of freedom of the lattice until it saturates and reaches infinite temperature. An alternative method of achieving this is through Floquet Heating. To show this we consider the closed Hubbard model under periodic driving in the form of a time-dependent inhomogeneous magnetic field with B(t) = V cos(Ωt) and f (i) describing the inhomogeneity. In the long-time limit the driving is expected to thermalise the system to infinite temperature [15,16]. However, because in this case the driving term commutes with all the generators of the η-symmetry we expect heating to occur only within the spin sector. In Fig. 3. we show how, by time-evolving an initial state under H(t), we realise long-time dynamics identical to that induced by the spin-dephasing. By resonantly driving the system, U = nΩ n ∈ Z + , the system thermalises quickly [17] and the observables are in strong agreement with both the grand-canonical description of (6) and long-time simulations of the map in Eq. (4). There is completely uniform off-diagonal long-range order in η-pairs (Fig. 3. Insets).
In Fig. 3a. we choose a magnetic field with a linear gradient f (i) = i, whilst in Fig. 3b. we choose a disordered field. In both cases the long-time dynamics are the same and we use this to emphasize that the choice of driving parameters and inhomogeneity is arbitrary. This choice only affects the time-scale on which the system relaxes and further data is presented in the Supplemental Material.
Conclusion -We have demonstrated that long-range ηpairing can be created and protected within the Hubbard model by directly heating the spinful degrees of freedom to infinite temperature. This destroys any coherences involving spinful particles which, in turn, creates uniform correlations between the η quasi-particles in the lattice with no identifiable length-scale.
Recent work has demonstrated that, with a judicious choice of parameters, applying a Gaussian pulse to the Mott-Insulating phase of the Hubbard model can transiently excite η-pairs [47]. This driving does not commute with the generators of the η-symmetry and hence completely uniform long-range correlations are not guaranteed. This does, however, open up the possibility of applying our heating scheme in conjuction with a sim-ilar optical driving process which can enhance the ηcorrelations in the system. Our heating would then continually spread these over all length-scales, creating a long-lived superconducting state.
We also anticipate further work using an alternate heating mechanism within the Hubbard model, or t-J model, which protects singlets and destroys other coherences. This could enhance superexchange pairing and induce superconductivity through nearest neighbour singlet pairing [12,48,49].
Finally, we emphasize the results of this work are due to the multiple SU(2) symmetries of the Hubbard model and not its microscopic details. Hence, ODLRO can be realised through dephasing in any model with 2 or more SU(2) symmetries (see Supplemental Material), as examples we suggest multi-band Hubbard models and the Richardson-Gaudin model [50] which also permit superconducting regimes.
We would like to thank S.R. Clark for our fruitful conversations and his beneficial input on optical modulation. We also thank C. The master equation we solve in the main text reads with For our numerical simulations we perform a stochastic unravelling of this master equation, known as the 'quantum trajectories' approach [36]. We use DMRG (Density Matrix Renormalisation Group) [37] to find the ground state of the Hamiltonian and then evolve this in time using the TEBD (Time Evolving Block Decimation) [38] algorithm. We focus on 1D lattice realisations of the model due to their numerical tractibility and the advantage provided by the aforementioned MPS (Matrix Product State) algorithms. We implemented these algorithms using the TNT (Tensor Network Theory Library) [39]. For all figures within the main text we use a bonddimension of χ = 1000 to ensure the corresponding SVD (Singular Value Decomposition) errors are minimal and the sum of the squares of the discarded singular values in a given time-step does not exceed = 1 × 10 −4 . Increasing the bond dimension has no effect on our results. For all simulations with a finite value of γ we perform N = 2000 trajectories to ensure convergence of the measured observables to within an uncertainty of 2%. We also checked that lowering the timestep δt from the value (δt = 0.01/τ ) used in our second-order Trotter decomposition of the propagator has no noticeable effect on our results.

USING SYMMETRY BASED DEPHASING TO INDUCE DISTANCE INVARIANT OFF-DIAGONAL CORRELATIONS
The result of engineering off-diagonal long-range correlations through dissipation is not specific to the Hubbard model nor the choice of Lindblad operators, it is a consequence of the multiple SU(2) symmetries of the Hamiltonian. To illuminate this we define a general Liouvillian map Λ in the manner of (8), with a Hamiltonian H and set of jump operators {L k }. We consider the situation when H has at least 2 SU(2) symmetries, i.e. there exists N ≥ 2 sets of operators {J ±,z i } i = 1, 2, ..., N where J ±,z i are the 3 generators for the SU(2) symmetry i of the Hamiltonian The commutativity between any pairs of generators from two different symmetries is encoded by the Kronecker delta. Any of the operators can be written as a sum of their projectors, via an eigendecomposition, i J − i due to their mutual eigenspaces. We now induce dephasing on the model in a symmetry protected manner by choosing {L k } such that for at least one symmetry j within i = 1, 2, ..., N . Provided that the map is unital ( k [L k , L † k ] = 0) then the commutation relations in Eq. (10) result in any projector associated with a generator satisfying (11) being a null eigenvector of the Liouvillian map Λ(P +−,z j,α ) = 0 [46]. The projectors P +−,z j,α for all the j satisfying (11), and those for the operators {C} relating to any remaining conserved quantities [51], fully span the kernel of Λ. It then follows that linear combinations (or any tracenormalised function which can be expanded as a power series) of these projectors is a steady state of the map Λ [46].
Let us write the steady states in this form where the sum is over all the projectors spanning the kernel of Λ. The following analysis is also valid if we write the steady state as any power-series expandable function of these projectors. Assuming we have a lattice model we define a permutation operator P i,i (P 2 i,i = 1) which exchanges sites i and i . If all the operators, {C}, whose projectors are involved in Eq. (12) are invariant under the same P i,i then Now let us choose one of the symmetries satisfying (11) and index it by k. The 3 generators for this symmetry are J +,−,z k and all of their projectors are contained in Eq. (12). We assume they can be written as a . Then we define the off-diagonal expectation value We now prove this is invariant under the choice of l and m by using the resolution of identity with the permutation operator P j,j as well as (14) Tr(ρ ss J + k,l J − k,m ) = Tr(ρ ss P 2 l,l J + k,l P 2 m,m J − k,m ) = Tr(ρ ss P l,l J + k,l P l,l P m,m J − k,m P m,m ) = Tr(ρ ss J + k,l J − k,m ) l = m, l = m , l = m , l = m.
We can easily choose the correct combination of indices on the permutation operators to lift the restrictions l = m and l = m and so we get: Hence, we show that these correlations have no identifiable length scale as they are uniform for any of the symmetries j, which satisfy (11), if we can find a valid two-site permutation operator which commutes with all the operators {O} = {J +−,z j } {C} at the same time. In the main text we choose spin dephasing L k = s z k to observe uniform, off-diagonal correlations in the η SU(2) symmetry of the Hubbard model. This dephasing means the η generators satisfy (11) but also that S z is an additional conserved quantity. Hence {O} = {η + , η − , η z , S z }, and the kernel of Λ is made up of the projectors of these operators (which are equivalent to the projectors of η + η − , N ↑ , N ↓ ). We can find a valid permutation operator for all of these operators simultaneously -it must swap two sites and carry over a change of sign on the local operators when it does this. We consider initial states (ground states of the Hamiltonian) which typically have decaying short-range correlations in η pairs. As η + η − is conserved these correlations will be spread over all length-scales of the system in order to satisfy Eq. (17). Hence, in our simulations, we see uniform offdiagonal pairing in the η sector.

STABILITY OF THE MASTER EQUATION UNDER PERTURBATION
We consider how the dynamics of our model is affected by perturbations in the dephasing. To simulate this we choose additional dephasing in the number basis of the system: as well as having {L k } = s z k with constant rate γ s , ∀k we add in another set of jump operators {L m } = n m,↑ + n m,↓ with constant rate γ c = 0.01γ s , ∀m.
Whilst in the long-time limit this additional dephasing completely destroys any coherences in the system, there exists a timescale in which the off-diagonal long-range order (ODLRO) appears before the effects of the additional dephasing become significant. Figure 4 shows that this additional dephasing does not prevent the observation of the ODLRO which forms (Fig. 4d) and is maintained prior to the system decaying to a thermal classical ensemble.
The window, indicated by the dashed lines in Figures  4a, 4b and 4c, shows the time it takes for the off-diagonal correlations to decay a quarter of their value at the first dashed line. The size of this window is not diminishing with system size. This suggests that even in the thermodynamic limit there is still a timeframe in which ODLRO can be observed when there is unwanted dephas- ing present. Additionally, numerics and perturbation arguments show that for small perturbations, the size of this window is controlled by the ratio γ s /γ c and so being able to tune γ s can help mitigate the effects of the unwanted dephasing. This is highlighted in the difference between Figs. 4b and 4d.

FLOQUET HEATING
In the main text we discussed that the properties of the steady states of (8) are reproduced by Floquet heating of the spin sector of the Hubbard model to infinite temperature. In order to show this we consider the Hubbard Hamiltonian with an additional time-dependent inhomogeneous magnetic field and f (i) describes the inhomogeneity of the field. In Figure. 5 we plot the convergence of the η-correlations to those expected from the Grand Canonical Ensemble (GCE), ρ GCE ∝ exp(µ 1 η + η − + µ 2 N ↑ + µ 3 N ↓ ), for 3 different choices of inhomogeneity. As expected the choice of f (i) is arbitrary and, along with the choice of parameters, only affects the transient dynamics and the rate of relaxation.

OSCILLATING COHERENCES
The Liouvillian in Eq. (8) also contains equally spaced imaginary eigenvalues [35] due to the existence of a ladder operator, η ± , of the Hamiltonian which also commutes with the Lindblad operators. The associated eigenvectors take the form of a raising of the steady state ρ ss : ρ nm ∝ (η + ) n ρ ss (η − ) m with the eigenvalues λ nm = 2i(m − n)µ. These states couple together sectors of the Hilbert Space with different particle numbers.
If one were able to initialise the system in a state with coherent superpositions between states with different particle numbers then, in the long-time limit, the dynamics of observables such as η x should be oscillatory [35]. In the main text we have focused solely on the properties of the steady state by enforcing a specific number of particles within our simulations.

QUASI-MOMENTUM DISTRIBUTIONS
We also plot the momentum distribution of the doublons following a quench from the ground state of the Hubbard model using the master equation (8). Figure 6. shows the evolution of the doublon structure factor D(qa, t) = 1 L L j,k=1 c † ↑,j c † ↓,j c ↓,k c ↑,k e i(k−j)qa (t), (19) in time for the full range of quasi-momenta q = 2πn/La n ∈ {0, 1, ..., L − 1}, with a the lattice spacing. In the long-time limit of the quench we see the flattening and equal excitation of all momenta modes, except for qa = π. This is characteristic of the steady states of the map in Eq. (8) due to their off-diagonal long-range order. The occupation of the momentum mode at qa = π is equal to η + η − /L, it is a constant of the evolution and for states with increasing values of η + η − the amplitude of this mode will grow (Fig. 6c), indicating the presence of a doublon condensate [52].