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

TL;DR

This paper demonstrates that heating the spin degrees of freedom in the Hubbard model to infinite temperature can induce long-range $ ext{eta}$-pairing correlations in the steady state, independent of the heating method, due to symmetry properties.

Contribution

It reveals how symmetry considerations enable long-range $ ext{eta}$-pairing in the Hubbard model's steady state under heating, extending the understanding of nonequilibrium states.

Findings

Steady states exhibit distance-invariant $ ext{eta}$-correlations.

Long-range $ ext{eta}$-pairing arises from symmetry, not microscopic details.

Heating to infinite temperature melts spin order but induces $ ext{eta}$-pairing.

Abstract

We show how, upon heating the spin degrees of freedom of the Hubbard model to infinite temperature, the symmetries of the system allow the creation of steady-states with long-range correlations between $\eta$-pairs. With induce this heating with either dissipation or periodic driving and evolve the system towards a nonequilibrium steady state, a process which melts all spin order in the system. The steady state is identical in both cases and displays distance-invariant off-diagonal $\eta$-correlations. These correlations were first recognised in the superconducting eigenstates in Yang's seminal paper [Phys. Rev. Lett 63, 2144 (1989)], which are a subset of our steady states. We show that our results are a consequence of symmetry properties and entirely independent of the microscopic details of the model and the heating mechanism.

Equations

$$H=-\tau \sum _{⟨ij⟩,\sigma }(c_{\sigma ,i}^{†}{c}_{\sigma ,j}+h.c)+U\sum _i{n}_{\uparrow ,i}{n}_{\downarrow ,i},$$

$${\eta }^{+}$$

$${\eta }^{-}$$

$${\eta }^z$$

$$\mathrm{Tr}(\rho {\eta }_i^{+}{\eta }_j^{-})=const.\forall i,ji\ne j,$$

$$\frac{\partial \rho }{\partial t}=\mathcal{L}\rho =-i[H,\rho ]+\gamma \sum _j$$

$$L_j$$

$$[H,A]=[L_j,A]=[L_j^{†},A]=0\forall j,$$

$${\rho }_{ss}\propto \exp ({\mu }_1{\eta }^{+}{\eta }^{-}+{\mu }_2{N}_{\uparrow }+{\mu }_3{N}_{\downarrow }),$$

$$H(t)=H+B(t)\sum _{i=1}^{L}f(i)s_i^z,$$

$$\frac{\partial \rho }{\partial t}=\mathcal{L}\rho =-i[H,\rho ]+\gamma \sum _{j=1}^M(L_j\rho L_j^{†}-\frac{1}{2}\{L_j^{†}{L}_j,\rho \}),$$

$$L_j=s_j^z=n_{\uparrow ,j}-n_{\downarrow ,j},$$

$$H=-\tau \sum _{⟨ij⟩,\sigma }(c_{\sigma ,i}^{†}{c}_{\sigma ,j}+h.c)+U\sum _{i=1}^M{n}_{\uparrow ,i}{n}_{\downarrow ,i}.$$

$$[H,J_i^z]$$

$$[J_i^z,J_{i^{\prime }}^{\pm }]$$

$$[L_k,J_j^{+-,z}]=0,$$

$${\rho }_{ss}\propto \sum _{\beta }{\lambda }_{\beta }{\mathcal{P}}_{\beta },$$

$$P_{i,i^{\prime }}O{P}_{i,i^{\prime }}=O⟹[O,P_{i,i^{\prime }}]=0⟹[{\mathcal{P}}_{\beta },P_{i,i^{\prime }}]=0\forall \beta .$$

$$[{\rho }_{ss},P_{i,i^{\prime }}]=0.$$

$$⟨J_{k,l}^{+}{J}_{k,m}^{-}⟩=\mathrm{Tr}({\rho }_{ss}{J}_{k,l}^{+}{J}_{k,m}^{-})l\ne m.$$

$$\mathrm{Tr}({\rho }_{ss}{J}_{k,l}^{+}{J}_{k,m}^{-})=\mathrm{Tr}({\rho }_{ss}{P}_{l,l^{\prime }}^2{J}_{k,l}^{+}{P}_{m,m^{\prime }}^2{J}_{k,m}^{-})=$$

$$\mathrm{Tr}({\rho }_{ss}{P}_{l,l^{\prime }}{J}_{k,l}^{+}{P}_{l,l^{\prime }}{P}_{m,m^{\prime }}{J}_{k,m}^{-}{P}_{m,m^{\prime }})=$$

$$\mathrm{Tr}({\rho }_{ss}{J}_{k,l^{\prime }}^{+}{J}_{k,m^{\prime }}^{-})l\ne m,l^{\prime }\ne m^{\prime },l\ne m^{\prime },l^{\prime }\ne m.$$

$$⟨J_{k,l}^{+}{J}_{k,m}^{-}⟩=⟨J_{k,l^{\prime }}^{+}{J}_{k,m^{\prime }}^{-}⟩l\ne m,l^{\prime }\ne m^{\prime }.$$

$$H(t)=H_{\mathrm{Hubbard}}(U,\tau )+B(t)\sum _{i=1}^{L}f(i)s_i^z,$$

$$B(t)=V\cos (\Omega t),$$

$$D(qa,t)=\frac{1}{L}\sum _{j,k=1}^L⟨c_{\uparrow ,j}^{†}{c}_{\downarrow ,j}^{†}{c}_{\downarrow ,k}{c}_{\uparrow ,k}{e}^{\mathrm{i}(k-j)qa}⟩(t),$$

$$H_L=H,$$

$$\displaystyle H_{I}=\int\big{(}\kappa_{\uparrow}\chi_{\uparrow}^{\dagger}(\mathbf{r})\chi_{\uparrow}(\mathbf{r})\phi^{\dagger}(\mathbf{r})\phi(\mathbf{r})+\ldots$$

$$\displaystyle\ldots\kappa_{\downarrow}\chi_{\downarrow}^{\dagger}(\mathbf{r})\chi_{\downarrow}(\mathbf{r})\phi^{\dagger}(\mathbf{r})\phi(\mathbf{r})\big{)}d\mathbf{r},$$

$$\displaystyle H_{\rm BEC}=\int\phi^{\dagger}(\mathbf{r})\bigg{[}-\frac{\nabla^{2}}{2m_{b}}+V_{\rm ext}(\mathbf{r})+\frac{g}{2}\phi^{\dagger}(\mathbf{r})\phi(\mathbf{r})\bigg{]}\phi(\mathbf{r})d\mathbf{r},$$

$${\partial }_t{\rho }_L(t)=-i[H_g,{\rho }_L(t)]+\dots$$

$$\displaystyle\ldots\sum_{l,l^{\prime}=1}^{M}f_{l,l^{\prime}}(t)\big{(}s^{z}_{l}s^{z}_{l^{\prime}}\rho_{L}(t)+\rho_{L}(t)s^{z}_{l}s^{z}_{l^{\prime}}-2s^{z}_{l}\rho_{L}(t)s^{z}_{l^{\prime}}\big{)},$$

$$H_g=H^{{}^{\prime }}-\sum _{l,l^{\prime }=1}^M{g}_{l,l^{\prime }}(t)s_l^z{s}_{l^{\prime }}^z,$$

Full text

Introduction - Driving and dissipation have recently emerged as transformative tools for dynamically evolving quantum systems into non-equilibrium phases with desirable properties Verstraete et al. [2009]; Cavalleri [2018]; Eisert et al. [2015]. In the context of strongly-correlated 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 Touzard et al. [2018]; Lüschen et al. [2017]. 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 Diehl et al. [2008]; Kraus et al. [2008]. 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 (d.o.f.) using terahertz laser pulses Hu et al. [2014]; Mitrano et al. [2016]. These optically driven systems are often modelled via Floquet driving - where the Hamiltonian is subject to a time-dependent periodic field Coulthard et al. [2017]; Sentef et al. [2016]; Knap et al. [2016]. 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 D’Alessio and Rigol [2014]; Lazarides et al. [2014]; Hermann et al. [2017] 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 that, counter-intuitively, the interplay between symmetry and heating can actually lead to the formation of steady states with coherent, long-range correlations. This heating can be achieved either with dephasing or with periodic driving and, as the generation of these correlations is based on symmetry, the microscopic details of the heating mechanism are unimportant. We believe this could mark an important step in explaining why such states are being observed in experiments Hu et al. [2014]; Mitrano et al. [2016] where the stringent conditions required in previous numerical work Coulthard et al. [2017]; Kaneko et al. [2019]; Sentef et al. [2016]; Knap et al. [2016] are not fulfilled.

In the aforementioned steady states we demonstrate how the formation of these long-range correlations results from the melting of all order in the complementary symmetry sector where the heating is applied. The competition between different types of order is a mechanism considered to underpin the formation of transient superconductivity Fausti et al. [2011]; Sentef et al. [2017] and thus we focus our work within this context.

Specifically, we apply driving or dephasing induced heating to the spin d.o.f. of bi-partite $D$-dimensional realisations of the Hubbard model, melting any order in this sector and reach robust mixed states with completely uniform long-range correlations between $\eta$-pairs. This pairing is known to provide a possible mechanism for superconductivity Yang [1962]; de Boer et al. [1995]; Essler et al. [1992]; Sentef et al. [2017]. Compared to previous proposals, which excite $\eta$-paired states through carefully tuned driving or dissipation Diehl et al. [2008]; Kaneko et al. [2019], our results are based on symmetry arguments. This means that the engineered steady states are independent of the model parameters and, for arbitrary initial states, guaranteed to have completely uniform $\eta$-correlations. In the case of dephasing we prove the results of our simulations by block-diagonalising the Liouvillian superoperator, giving an explicit parametric form for the steady states.

Model - The Hubbard model, in second quantised form, reads

[TABLE]

where $c_{\sigma,i}^{\dagger}$ and its adjoint are the usual creation and annihilation operators for a fermion of spin $\sigma\in\{\uparrow,\downarrow\}$ on site $i$. Additionally, $n_{\sigma,i}$ is the number operator for a particle of spin $\sigma$ on site $i$ and $\tau$ and $U$ play the role of kinetic and interaction energy scales respectively. In this work we set $\hbar=1$.

The Hamiltonian in Eq. (1) has a rich structure due to the two SU(2) symmetries it possesses Essler et al. [2005]. The first of these is the ‘spin’ symmetry which relates to spinful particles (singlons) $\sigma\in\{\uparrow,\downarrow\}$. The second, often referred to as ‘$\eta$-symmetry’, is central to this letter and relates to spinless quasi-particles (doublons and holons) $\sigma\in\{\uparrow\downarrow,{\rm vac}\}$. It can be interpreted as a type of particle-hole symmetry. This $\eta$-symmetry is revealed by introducing the associated operators

[TABLE]

with $\eta^{+}_{i}$ ($\eta^{-}_{i}$) describing the creation (annihilation) of a doublon on site $i$ with an alternating site-dependent phase. The operators in Eq. (2) fulfill the relations $[H,\eta^{\pm}]=0$ and $[H,\eta^{+}\eta^{-}]=[H,\eta^{z}]=0$ and commute with all the generators of the spin symmetry.

The presence of $\eta$-pairing superconductivity in the Hubbard model is a phenomenon first recognised by C. N. Yang in his seminal paper Yang [1989]. There it was proved that the pure states $\ket{\psi}\propto(\eta^{+})^{N}\ket{{\rm vac}}$ are eigenstates of $H$ and possess off-diagonal long-range order (ODLRO) in the form of doublon-doublon correlations

[TABLE]

where $\rho=\ket{\psi}\bra{\psi}$. 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 Yang [1962]; Sewell [1990]; Nieh et al. [1995]. 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. By driving the Hubbard model in the spin basis these are destroyed and we can consistently engineer states with long-range uniform correlations in $\langle\eta^{+}_{i}\eta^{-}_{j}\rangle$.

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 Kruchinin et al. [2010]; Monthoux and Pines [1992]; Monthoux et al. [1991]. Our dephasing mimics these scattering events.

Despite the ‘toy’ nature of our model the Hubbard Hamiltonian can be accurately realised by loading ultracold fermionic atoms into optical lattices Schneider et al. [2012]; Jaksch and Zoller [2005]; Lewenstein et al. [2007]. These quantum simulators offer precise experimental control over the microscopic details of the system. In the Supplemental Material (SM) Note [1] we show how dephasing can occur by immersing the lattice into a homogeneous Bose-Einstein condensate Klein et al. [2007]; Bruderer et al. [2007]. If the interactions are tuned with Feshbach resonances Theis et al. [2004]; Inouye et al. [1998] so that the scattering amplitudes between the two fermionic spin states and the bosons are equal in magnitude and opposite in sign, then dephasing will occur solely in the spin sector.

Results - We couple the Hubbard Hamiltonian in Eq. (1) 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

[TABLE]

with Lindblad operators $s^{z}_{j}$ on each site $j$. Because these operators are restricted to the spin SU(2) symmetry sector we obtain $[L_{j},\eta^{\pm}]=[L_{j},\eta^{z}]=0\ \forall j$.

For Eq. (4) any operator $A$ which satisfies

[TABLE]

is a null eigenvector of $\mathcal{L}$ and the adjoint map $\mathcal{L}^{\dagger}$, which means the associated observable $\langle A\rangle$ is conserved. Thus, the map in Eq. (4) conserves $\langle\eta^{+}\eta^{-}\rangle$, $\langle\eta^{z}\rangle$ and $\langle S^{z}\rangle=\sum_{i}\langle s^{z}_{i}\rangle$. Through a power series expansion, any appropriately trace normalised function $f(\eta^{+}\eta^{-},\eta^{z},S^{z})$ of these operators is a steady state of $\mathcal{L}$, i.e. if $\rho_{ss}=f(\eta^{+}\eta^{-},\eta^{z},S^{z})$ then $\mathcal{L}\rho_{ss}=0$. In the SM Note [1] we show that steady states of this form have ${\rm Tr}(\rho_{ss}\eta^{+}_{i}\eta^{-}_{j})={\rm const}\ \forall i,j$ (see also Buca et al. [2019]). We also show how this result of heating-induced long-range order can be extended to arbitrary models with multiple SU(2) symmetries. We emphasize that the expected steady state $\rho_{ss}$ is realisable for any bi-partite $D$ dimensional realisation of the spin-dephased Hubbard model.

Here, we focus on 1D lattice realisations due to their numerical tractability. We perform calculations where the system is initialised in a given state and then evolved by solving the time-dependent master equation in Eq. (4) and quenching $U$, providing a direct comparison with the dynamics of the closed system. In the open system the value of $U$ quenched to solely sets the time-scale of relaxation and long-range order will always appear. We use a quantum trajectories approach Daley [2014] to solving Eq. (4) combined with Density Matrix Renormalisation Group White [1992] and Time Evolving Block Decimation Vidal [2003] algorithms for finding the ground state and evolving the wavefunction of the system respectively. These simulations were performed with the Tensor Network Theory library Al-Assam et al. [2017]. Further numerical details and analysis can be found in the SM Note [1].

In Fig. 1. we plot, for various distances $j$, the averaged quantity $\braket{\eta^{+}_{i}\eta^{-}_{i+j}}$ over time for a quench from the ground state of the Hubbard Hamiltonian. We also include the matrix of correlations of $\langle\eta^{+}_{i}\eta^{-}_{j}\rangle$ for the density matrix following the quench. For comparison, we show the case where the system is not coupled to the environment, $\gamma=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 $\langle\eta^{+}_{i}\eta^{-}_{j}\rangle$ is finite and constant for all $i\neq j$. We observe that this is facilitated by a decrease in the short-range $\eta$-pairing correlations in order to allow the long-range ones to increase: the non-equilibrium dynamics involve a ‘spreading’ of the $\eta$-correlations over all length scales of the system which is necessary due to the conservation of $\eta^{+}\eta^{-}$. Plots of the corresponding doublon momentum distribution are discussed in the SM Note [1].

We also show 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, breaking the strong symmetry relations and causing the steady state of the Liouvillian to contain no coherences. Nonetheless, we set the strength of the charge-dephasing to values around $1\%$ of the spin-dephasing and find there exists an intermediate time-scale where the results of Fig. 1 can be observed. We also find that 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: the dissipative evolution induced by the charge-dephasing becomes increasingly frozen out as this ratio increases. In particular, our data suggests that the length of this window can be extended by increasing the coupling of the spin dephasing whilst the other parameters in the model remain constant.

To obtain a more intuitive idea of the nature of the steady state $\rho_{ss}=f(\eta^{+}\eta^{-},\eta^{z},S^{z})$ we consider a physically motivated parametrisation of the steady state function

[TABLE]

where we have exchanged $\eta^{z}$ and $S^{z}$ with $N_{\uparrow}=\sum_{i}n_{\uparrow,i}$ and $N_{\downarrow}=\sum_{i}n_{\downarrow,i}$ as they can be expressed as linear combinations of these operators. Calculations on small lattices show that the parametrisation in Eq. (6) captures all relevant observables when compared to numerical simulations which reach the long-time limit in Eq. (4). Equation (6) describes a generalized grand-canonical-like equilibrium state (GCE) with the Lagrange multipliers $\mu_{1},\mu_{2},\mu_{3}$ associated to each of the conserved quantities. Notably, the Lagrange multiplier for the Hamiltonian is [math], i.e. at infinite temperature. As is typical of states in this form the Lagrange multipliers are independent of the quench parameters and are governed solely by the values of the conserved quantities $\langle\eta^{+}\eta^{-}\rangle$, $\langle N_{\uparrow}\rangle$ and $\langle N_{\downarrow}\rangle$ associated with the initial state. The existence of a ‘quadratic’ charge ($\eta^{+}\eta^{-}$) in the GCE state is crucial for the presence of non-trivial long-range order and sets it apart from Generalised Gibbs Ensemble states Ilievski et al. [2015]; Rigol et al. [2007]; Vidmar and Rigol [2016]; Mur-Petit et al. [2018].

To emphasize the properties of the stationary state in Eq. (6) we plot, in Fig. 2., the projection of the density matrix in the spin ($s$) and $\eta$ ($d$) d.o.f. $\mathcal{P}_{s,d}\rho\mathcal{P}_{s,d}$. Explicitly, these projectors remove any basis vectors which contain, respectively, doublons ($\sigma\in\{\uparrow\downarrow,\rm vac\}$) and singlons ($\sigma\in\{\uparrow,\downarrow\}$) 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 Buča and Prosen [2012] $\langle\eta^{+}\eta^{-}\rangle$, $\langle N_{\uparrow}\rangle$ and $\langle N_{\downarrow}\rangle$ can be used to block-diagonalise the Liouvillian and hence the degeneracy of $\rho_{ss}$ is determined by the distinct combinations of the eigenvalues of these operators. The steady states are the states with maximum $\eta$-order in each block, quantified by the amplitude of the uniform off-diagonal correlations $|\langle\eta^{+}_{i}\eta^{-}_{j}\rangle|,\ \ i\neq j$. Notably, Yang’s states $\ket{\psi}\propto(\eta^{+})^{N}\ket{{\rm vac}}$ form a subset of these steady states. They are the only pure steady states as they exist in the decoherence free subspace Lidar et al. [1998] of Eq. (4).

The infinite temperature ensemble pictured in Fig. 2c. shows that the dephasing has continuously pumped energy into the spin d.o.f. 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

[TABLE]

with $B(t)=V\cos(\Omega t)$ and $f(i)$ describing the inhomogeneity. In the long-time limit the driving is expected to thermalise the system to infinite temperature D’Alessio and Rigol [2014]; Lazarides et al. [2014]. However, because in this case the driving term commutes with all the generators of the $\eta$-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\Omega\ n\in\mathbb{Z}^{+}$, the system thermalises quickly Hermann et al. [2017] and the observables are in good 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 $\eta$-pairs (Fig. 3. Insets).

In Fig. 3a. we choose a linear magnetic field $f(i)=i$, whilst in Fig. 3b. we choose a disordered field. In both cases the long-time dynamics are the same and this emphasizes that the choice of driving parameters and inhomogeneity is arbitrary - they only affect the time-scale on which the system relaxes.

Conclusion - We have demonstrated that long-range $\eta$-pairing can be created and protected within the Hubbard model by directly heating the spinful d.o.f. to infinite temperature. This destroys any coherences involving spinful particles which, in turn, creates uniform long-range correlations between the $\eta$ quasi-particles in the lattice.

Recent work has demonstrated that, with a judicious choice of parameters, applying a Gaussian pulse of a specific duration to the Mott-Insulating phase of the Hubbard model can excite $\eta$-pairs Kaneko et al. [2019]. This driving does not commute with the generators of the $\eta$-symmetry and hence uniform long-range correlations are not guaranteed. This does, however, open up the possibility of applying our heating scheme in conjuction with a similar driving process which can enhance the $\eta$-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 Coulthard et al. [2017]; Hirsch [1985]; Kotliar and Liu [1988].

Finally, we emphasize the results of this work are due to the multiple ${\rm SU}(2)$ symmetries of the Hubbard model and not its microscopic details. Hence, ODLRO can be realised through dephasing in any model with multiple symmetries (see SM Note [1]), as examples we suggest multi-band Hubbard models and the Richardson-Gaudin model Dukelsky et al. [2004] which also permit superconducting regimes.

I Supplementary Material To Heating-Induced Long-Range $\eta$-Pairing in the Hubbard Model

II Numerical Procedure

The master equation we solve in the main text reads

[TABLE]

with

[TABLE]

For our numerical simulations we perform a stochastic unravelling of this master equation, known as the ‘quantum trajectories’ approach Daley [2014]. We use DMRG (Density Matrix Renormalisation Group) White [1992] to find the ground state of the Hamiltonian and then evolve this in time using the TEBD (Time Evolving Block Decimation) Vidal [2003] algorithm. We focus on 1D lattice realisations of the model due to their numerical tractibility. We implemented these algorithms using the TNT (Tensor Network Theory) Library Al-Assam et al. [2017]. For all figures within the main text we use a bond-dimension of $\chi=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 $\epsilon=1\times 10^{-4}$. Increasing the bond dimension has no effect on our results. For all simulations with a finite value of $\gamma$ 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 $\delta t$ from the value ($\delta t=0.01/\tau$) used in our second-order Trotter decomposition of the propagator has no noticeable effect on our results.

Despite using a TEBD approach we are not able to reach particularly large system sizes as the states reached in the long-time limit have both completely long-range entanglement and a significant amount of classical correlations. This can cause the Matrix Product State representation to be highly inefficient Datta and Vidal [2007] and we find the bond dimensions required for our simulations, at least at the levels of the individual trajectories, is large and prevents us from being able to make effective use of the aforementioned Matrix Product algorithms.

In Fig. 4 we demonstrate this: the cumulative total of the squares of the discarded singular values (a measure of the accuracy of the simulation) over time for an average trajectory in the open system is several orders of magnitude larger than that of the corresponding closed system. Unless a very large bond dimension is used this results in a critical failure of the numerics to conserve physical symmetries of the system, such as $\langle\eta^{+}\eta^{-}\rangle$ - which should, for each trajectory, always remain [math] in the example considered.

III 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 $\Lambda$ in the manner of (II), 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\geq 2$ sets of operators $\{J_{i}^{\pm,z}\}\quad i=1,2,...,N$ where $J_{i}^{\pm,z}$ are the $3$ generators for the SU(2) symmetry $i$ of the Hamiltonian

[TABLE]

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,

$J_{i}^{+-,z}=\sum_{\alpha=1}^{D_{i}^{+-,z}}\lambda_{i,\alpha}^{+-,z}\mathcal{P}_{i,\alpha}^{+-,z}$. We have combined $J_{i}^{\pm}$ to form $J_{i}^{+-}=J_{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

[TABLE]

for at least one symmetry $j$ within $i=1,2,...,N$. Provided that the map is unital ($\sum_{k}[L_{k},L^{\dagger}_{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 $\Lambda(\mathcal{P}_{j,\alpha}^{+-,z})=0$ Buča and Prosen [2012]. The projectors $\mathcal{P}_{j,\alpha}^{+-,z}$ for all the $j$ satisfying (11), and those for the operators $\{C\}$ relating to any remaining conserved quantities Albert et al. [2016], fully span the kernel of $\Lambda$. It then follows that linear combinations (or any trace-normalised function which can be expanded as a power series) of these projectors is a steady state of the map $\Lambda$ Buča and Prosen [2012].

Let us write the steady states in this form

[TABLE]

where the sum is over all the projectors spanning the kernel of $\Lambda$. 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^{\prime}}$ ($P_{i,i^{\prime}}^{2}=1$) which exchanges sites $i$ and $i^{\prime}$. If all the operators, $\{O\}=\{J_{j}^{+-,z}\}\ \bigcup\ \{C\}$, whose projectors are involved in Eq. (12) are invariant under the same $P_{i,i^{\prime}}$ then

[TABLE]

It immediately follows that

[TABLE]

Now let us choose one of the symmetries satisfying (11) and index it by $k$. The $3$ generators for this symmetry are $J_{k}^{+,-,z}$ and all of their projectors are contained in Eq. (12). We assume they can be written as a sum of purely local operators on the $M$ site lattice: $J_{k}^{+,-,z}=\sum_{l=1}^{M}J_{k,l}^{+,-,z}$. Then we define the off-diagonal expectation value

[TABLE]

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^{\prime}}$ as well as (14)

[TABLE]

We can easily choose the correct combination of indices on the permutation operators to lift the restrictions $l\neq m^{\prime}$ and $l^{\prime}\neq m$ and so we get:

[TABLE]

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_{j}^{+-,z}\}\ \bigcup\ \{C\}$ at the same time.

In the main text we choose spin dephasing $L_{k}=s_{k}^{z}$ to observe uniform, off-diagonal correlations in the $\eta$ SU(2) symmetry of the Hubbard model. This dephasing means the $\eta$ generators satisfy (11) but also that $S^{z}$ is an additional conserved quantity. Hence $\{O\}=\{\eta^{+},\eta^{-},\eta^{z},S^{z}\}$, and the kernel of $\Lambda$ is made up of the projectors of these operators (which are equivalent to the projectors of $\eta^{+}\eta^{-}$, $N_{\uparrow}$, $N_{\downarrow}$). 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 (e.g. ground states of the Hamiltonian) which typically have decaying short-range correlations in $\eta$ pairs. As $\langle\eta^{+}\eta^{-}\rangle$ 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 off-diagonal correlations in the $\eta$ sector.

IV 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, ‘unwanted’, dephasing in the number basis of the system - as well as having $\{L_{j}\}=s_{j}^{z}$ with constant rate $\gamma^{s}\ \ \forall j$ (see Eq. (II)) we add in another set of local jump operators $\{L_{m}\}=n_{m,\uparrow}+n_{m,\downarrow}$ with constant rate $\gamma^{c}=\lambda\gamma^{s}\ \ \forall m$. The indices $j$ and $m$ each run over the lattice sites in the system. We have used local number dephasing because this additional dephasing breaks the strong symmetry relation $[\eta^{+}\eta^{-},L_{m}]\neq 0\ \forall m$. Consequently, in the long-time limit this additional dephasing completely destroys any coherences in the system and the steady state is a featureless infinite temperature ensemble.

We set the amplitude of the unwanted dephasing to be around $2$ orders of magnitude less than the spin-dephasing. In cold atomic systems parameters can be tuned widely and unwanted decoherence mechanisms strongly suppressed Lewenstein et al. [2013]. In quantum materials microscopic models are not as well established. However, the rates of different decoherence mechanisms can vary widely and careful engineering of these materials can significantly reduce their sensitivity to specific sources of decoherence Awschalom et al. [2013]. As a result decoherence rates can easily differ by orders of magnitude, in line with the assumptions in our work.

Figure 5 shows that, for $\lambda=\gamma^{c}/\gamma^{s}=0.01$, the additional dephasing does not prevent the observation of the ODLRO which forms and is maintained prior to the system decaying to a thermal classical ensemble. The window, indicated by the dashed lines, shows the time it takes for the off-diagonal correlations to decay to $1/3$ of their value at the first dashed line. The data in Fig. 5 suggests the size of this window is not diminishing with system size. This is consistent with previous research studying the rate of relaxation in bulk-dephased many-body systems Žnidarič [2015] or in random Liouvillians Can [2019] - they find this rate either decays or saturates with system size.

Additionally, Fig. 6 shows that for small perturbations this window can be lengthened by increasing the ratio $\gamma^{s}/\gamma^{c}$. Thus, in an experimental setup, being able to tune $\gamma_{s}$ could help mitigate the effects of the unwanted dephasing. Physically, the increased amplitude of the wanted (spin) dephasing compared to the unwanted (number) dephasing means the system is ‘measured’ with increasing frequency by the spin dephasing in comparison to the frequency with which it is ‘measured’ by the number dephasing. As the ratio of these amplitudes increases then the spin dephasing effectively arrests the time evolution induced by the number-dephasing. Repeated measurements in the spin basis are freezing the system’s dynamics with respect to the dissipative evolution caused by the unwanted dephasing.

V Floquet Heating

In the main text we discussed that the properties of the steady states of (II) 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

[TABLE]

and $f(i)$ describes the inhomogeneity of the field. In Figure. 7 we plot the convergence of the $\eta$-correlations to those expected from the Grand Canonical Ensemble (GCE), $\rho_{\rm GCE}\propto\exp(\mu_{1}\eta^{+}\eta^{-}+\mu_{2}N_{\uparrow}+\mu_{3}N_{\downarrow})$, for $3$ different choices of inhomogeneity. As expected the choice of $f(i)$ is arbitrary and, along with the choice of system parameters, only affects the transient dynamics and the rate of relaxation.

VI Oscillating Coherences

The Liouvillian in Eq. (II) also contains equally spaced imaginary eigenvalues Buca et al. [2019] due to the existence of a ladder operator, $\eta^{\pm}$, of the Hamiltonian which also commutes with the Lindblad operators. The associated eigenvectors take the form of a raising of the steady state $\rho_{ss}$: $\rho_{nm}\propto(\eta^{+})^{n}\rho_{ss}(\eta^{-})^{m}$ with the eigenvalues $\lambda_{nm}=2{\rm i}(m-n)\mu$. 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 $\eta^{x}$ should be oscillatory Buca et al. [2019]. 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.

VII 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 (II). Figure 8. shows the evolution of the doublon structure factor

[TABLE]

in time for the full range of quasi-momenta $q=2\pi n/La$  $n\in\{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=\pi$. This is characteristic of the steady states of the map in Eq. (II) due to their off-diagonal long-range order. The occupation of the momentum mode at $qa=\pi$ is equal to $\langle\eta^{+}\eta^{-}\rangle/L$, it is a constant of the evolution and for states with increasing values of $\langle\eta^{+}\eta^{-}\rangle$ the amplitude of this mode will grow (Fig. 8c), indicating the presence of a doublon condensate Penrose and Onsager [1956].

VIII Implementing Spin Dephasing in a Cold Atom Setup

We now discuss the possibility of simulating the spin-dephased Hubbard model in a cold-atom setup. The Hamiltonian in Eq. 9 can be accurately realised by loading an ultracold gas of fermionic atoms into a lattice structure created with counter-propagating laser beams Jaksch and Zoller [2005]; Lewenstein et al. [2007, 2013]. These optical lattices provide the experimentalist precise control over the microscopic parameters of the system.

In order to engineer dephasing solely in the spin degrees of freedom of the lattice we consider the possibility of immersing the lattice into a homogeneous Bose-Einstein Consendate (BEC), which was also discussed in Buca et al. [2019] in the context of achieving dephasing in the number basis of the Hubbard model. The interactions between the BEC and the lattice atoms create deformations in the condensate which can be described as polarons, or coherent states of phonons Klein et al. [2007].

To explicitly derive the dynamics of the lattice in this environment we closely follow methodology of Klein et al. [2007] and Bruderer et al. [2007]. We start by writing the Hamiltonian $H_{\rm tot}=H_{L}+H_{I}+H_{\rm BEC}$ of the total system

[TABLE]

with $H_{L},H_{I}$ and $H_{\rm BEC}$ being the lattice (see Eq. 9), interaction and BEC Hamiltonians respectively. The operators $\chi_{\uparrow}$ and $\chi_{\downarrow}$ are the field operators of the two fermionic levels which couple to the bosonic field ($\phi^{\dagger}(\mathbf{r})$) with amplitudes $\kappa_{\uparrow}$ and $\kappa_{\downarrow}$ respectively. The mass of a condensate atom is given by $m_{b}$, $V_{\rm ext}$ is an external trapping potential for the BEC and $g$ is the interaction strength between the condensate atoms.

In order for the condensate to affect the dynamics of the lattice solely in the spin sector we set the coupling strengths to be equal and opposite, i.e. $\kappa_{\uparrow}=-\kappa_{\downarrow}=\kappa$. This could be achieved via Feshbach resonances Theis et al. [2004]; Inouye et al. [1998] which can be used to tune atomic scattering lengths over a wide range of values. It is also necessary that $\kappa$ is sufficiently small $\kappa\ll gn^{0}\epsilon^{D}$, where $n_{0}$ is the BEC density in the centre of the trap, $\epsilon$ is the healing length and $D$ is the system dimension Klein et al. [2007]. By solving the Gross-Pitaevskii equation and treating the resulting deformations in the BEC as coherent states of phonons an effective, discretized, Hamiltonian for the system can be derived. We then trace out the phonon degrees of freedom, invoking the rotating wave and Born approximations, to derive a Master equation for the density matrix of the lattice $\rho_{L}(t)$

[TABLE]

where $H^{{}^{\prime}}$ is the Hubbard Hamiltonian (the parameters have been slightly modified due to the BEC presence) and $f_{l,l^{\prime}}(t)$ and $g_{l,l^{\prime}}(t)$ are a pair of time-dependent, short-range (based on ‘typical’ parameters for a BEC) functions which become time-independent in the long-time limit. Notably, however, the explicit form of these functions does not alter the strong symmetries of the system (they are still $\eta^{+}\eta^{-}$, $S^{z}$ and $N_{\rm total}$). Therefore the steady states of this Master equation are guaranteed to have long-range $\eta$-paired correlations in the long-time limit.