Dynamical preparation of a steady ODLRO state in the Hubbard model with local non-Hermitian impurity

In recent years, nonequilibrium dynamics induced by dissipation exhibits intriguing properties. As an effective description, the non-Hermitian Hamiltonian arises when the system experiences dissipation to an environment Daley (2014); Dalibard et al. (1992). Recent years have seen a growing interest in non-Hermitian descriptions of condensed-matter systems which have not only extended the domain of condensed-matter physics with inspiring insights Lee (2016); Kunst et al. (2018); Yao et al. (2018); Gong et al. (2018); El-Ganainy et al. (2018); Nakagawa et al. (2018); Shen and Fu (2018); Wu et al. (2019); Yamamoto et al. (2019); Song et al. (2019); Yang and Hu (2019); Hamazaki et al. (2019); Kawabata et al. (2019a, b); Lee et al. (2019); Yokomizo and Murakami (2019); Jin et al. (2020) but also provided a fruitful framework to elucidate inelastic collisions Xu et al. (2017), disorder effects Shen and Fu (2018); Hamazaki et al. (2019), and system-environment couplings Nakagawa et al. (2018); Yang and Hu (2019); Song et al. (2019). In particular, the interplay between non-Hermiticity and interaction can give rise to exotic quantum many-body effect, ranging from non-Hermitian extensions of Kondo effect Nakagawa et al. (2018); Lourenço et al. (2018), many-body localization Hamazaki et al. (2019), Fermi surface in coordinate space Mu et al. (2020), to fermionic superfluidity Yamamoto et al. (2019); Okuma and Sato (2019). It has been shown that the cooperation between the non-Hermiticiy and interaction can alter drastically the macroscopic behavior that has been established in the Hermitian physics. Meanwhile, the generation of the superconducting-like states that have been induced by dissipation has received great attention Diehl et al. (2008); Kraus et al. (2008); Sentef et al. (2016); Mitrano et al. (2016); Coulthard et al. (2017). However, these schemes rely on judicious engineering of the system parameters to avoid thermalization.

Exceptional points (EPs) are degeneracies of non-Hermitian operators where the corresponding eigenstates coalesce into one state leading to a incomplete Hilbert space Berry (2004); Heiss (2012); Miri and Alù (2019); Zhang and Gong (2020).The peculiar features around EP have sparked tremendous attention to the classical and quantum photonic systems. The corresponding intriguing dynamical phenomena include asymmetric mode switching Doppler et al. (2016), topological energy transfer Xu et al. (2016), robust wireless power transfer Assawaworrarit et al. (2017), and enhanced sensitivity Wiersig (2014, 2016); Hodaei et al. (2017); Chen et al. (2017) depending on their EP degeneracies. Notably, the high-order EP with $N$ coalescent states (EPN) attract much more interest recently. Many works have been devoted to the formation of the EPN and corresponding topological characterization in both theoretical and experimental aspects Ding et al. (2016); Xiao et al. (2019); Pan et al. (2019); Zhang et al. (2020).

Given the above two rapidly growing fields, namely, the non-Hermitian interacting system and dynamics of the EP, we are motivated to examine how non-Hermiticity impacts on the interacting system, especially the dynamics when the high-order EPs of the interacting system presents. In this paper, we uncover the effect of cooperation between particle dissipation and pair fluctuation on the strongly-correlated system by concentrating on the non-Hermitian Hubbard model. We show that the considered interacting system can respect the full real spectrum even no obvious symmetry presents. Remarkably, we find that a balance local non-Hermitian impurity can induce the formation of high-order EPs in the spectrum in the way that degenerate states with different symmetry of parent Hermitian system coalesce. The $\eta$-pairing mechanism plays a vital role to achieve this intriguing property. Based on the performance of the system at EP, a scheme that produces a nonequilibrium steady superconducting-like state is proposed. Specifically, for an arbitrary initial state even a vacuum state, the critical system can drive it to the coalescent state that favors superconductivity manifested by the off-diagonal long-range order (ODLRO). Such a dynamical scheme can be realized no matter where the local coupling is applied and the correlation of the final steady state is independent of the relative distance between the two sites. Therefore, our finding is distinct from the previous investigations Tindall et al. (2019); Kaneko et al. (2019), and offer an alternative mechanism for generating superconductivity through nonequilibrium dynamics. On the other hand, the remarkable observation from our work can trigger further studies of both fundamental aspects and potential applications of critical non-Hermitian strongly correlated systems.

The rest of this paper is organized as followed. In Section  II, we present the general property of the considered non-Hermitian Hamiltonian especially focuses on the mechanism of the formation of high-order EP. Section  III is devoted to demonstrating how a local coupling can make the degenerate eigenstate coalesce. In Section  IV, the generation of the superconducting state based on the critical non-Hermitian Hubbard model is proposed. Section  V concludes this paper. Some nonessential details of our calculation are placed in the Appendix.

We consider a non-Hermitian extension of Hubbard model for two-component fermions. The Hamiltonian is in the form

where $\bm{h}_{j}=g_{j}\left(\lambda,\text{ }0,\text{ }i\gamma\right)$ with $\left\{{g_{j}}\right\}$ describing inhomogeneity and a set of arbitrary numbers, and the corresponding $\eta$-operators are defined by

obeying the Lie algebra, i.e., $[\eta_{j}^{+},$ $\eta_{j}^{-}]=2\eta_{j}^{z}$ and $[\eta_{j}^{z},$ $\eta_{j}^{\pm}]=\pm\eta_{j}^{\pm}$. Here $c_{i,\sigma}$ is the annihilation (creation) operator for an electron at site $i$ with spin $\sigma$ and $n_{i,\sigma}=c_{i,\sigma}^{\dagger}c_{i,\sigma}$. $H_{0}$ is a standard Hubbard model on a bipartite lattice where $t_{ij}$ and $U$ play the role of kinetic and interaction energy scale; $H_{I}$ describes the non-Hermitian impurity that consists of on-site pair fluctuation and imaginary magnetic field, which can be achieved by coupling the system to the environment and are within the reach of both the ultracold atom Lee and Chan (2014) and photonic experiments Fausti et al. (2011); Hu et al. (2014); Kaiser et al. (2014); Mitrano et al. (2016); Cantaluppi et al. (2018).

Now we turn to investigate the symmetry of the considered model. The Hamiltonian in Eq. (2) has a rich structure due to the two sets of SU(2) symmetries it possesses. The first, often referred to as $\eta$ symmetry, is central to this paper and can be characterized by the generators $\eta^{\pm}=\sum_{j}\eta_{j}^{\pm}$ and $\eta^{z}=\sum_{j}\eta_{j}^{z}$. It relates to spinless quasiparticles (doublons and holons) and can be interpreted as a type of particle-hole symmetry. The second of these is spin symmetry, and the corresponding generator $\bm{s}$ can be obtained by replacing $c_{j,\downarrow}$ of Eqs. (4)-(6) with $\left(-1\right)^{j}c_{j,\downarrow}^{\dagger}$, that is

It can be readily seen that the operators in Eqs. (4)-(6) fulfill the relations $\left[H,\eta^{\pm}\right]=\left[H,\eta^{z}\right]=0$ and commute with all the generators of the spin symmetry. The presence of $H_{I}$ spoils two such SU(2) symmetries. However, the two Hamiltonians $H_{0}$ and $H_{I}$ commute with each other when the interacting strength $g_{j}$ is homogeneous such that $H_{I}$ can be treated as $g\left(\lambda\eta^{x}+i\gamma\eta^{z}\right)$. As such the two Hamiltonians shares the common eigenstates. Although the two Hamiltonians commutes with each other, $H_{I}$ profoundly changes the energy spectrum of the whole system due to the emergence of EPN, which is distinct from the Hermitian case. In this paper, we focus on the subspace spanned by the $\eta$-pairing state $\left|\psi_{N_{\eta}}\right\rangle=\Omega^{-1}\left(\eta^{+}\right)^{N_{\eta}}\left|\mathrm{Vac}\right\rangle$, where $\left|\mathrm{Vac}\right\rangle$ is a vacuum state with no electrons and $N_{\eta}$ is the number of $\eta$ pairs. What makes the state $\left|\psi_{N_{\eta}}\right\rangle$ is special is the fact that it has been shown to have off-diagonal long-range order (ODLRO) in the form of doublon-doublon correlations, $\langle\psi_{N_{\eta}}|\eta_{i}^{\dagger}\eta_{j}^{-}\left|\psi_{N_{\eta}}\right\rangle=\mathrm{const}$, ($i\neq j$). The nonzero value of such quantity implies both the Meissner effect and flux quantization and therefore provides a possible definition of superconductivity Yang (1962). These states are the $2N+1$ fold degenerate eigenstates of both $H_{0}$ and $\eta^{2}$. Yang (1989) Hence, all these states can be expressed as $\{\left|N,l\right\rangle\}$, where $N$ and $l$ are associated with the eigenvalues of $\eta^{2}$ and $\eta^{z}$, i.e.,

with $\left|N,l\right\rangle=\Omega^{-1}\left(\eta^{+}\right)^{N+l}\left|\mathrm{Vac}\right\rangle$ and $\Omega=\sqrt{C_{2N}^{N+l}}$. For the homogeneous $g_{j}=g$, the matrix of $H$ in the doublon invariant subspace has the form

with $J_{l}=\sqrt{(N-l)(N+l+1)}$. It describes a $\mathcal{PT}$ -symmetric hypercube graph of $2N+1$ dimension Zhang et al. (2012). The EP at $\left|\lambda\right|=\left|\gamma\right|$ divides the system into two different phases: $\mathcal{PT}$ -symmetric broken ($\left|\lambda\right|<\left|\gamma\right|$) and unbroken regions ($\left|\lambda\right|>\left|\gamma\right|$). The unique feature of $H_{\mathrm{doub}}$ is that all the eigenstates coalesce at $\left|\lambda\right|=\left|\gamma\right|$ forming a high-order EP with order of $2N+1$. The corresponding coalescent state is

with Dirac normalization. If the system is initialized in this invariant subspace, then all the dynamical property of the system is solely determined by the non-Hermitian Hamiltonian $H_{\mathrm{doub}}$. Specifically, when $H_{\mathrm{doub}}$ respects full real spectrum, the physical observable will show a periodic oscillating behavior. This dynamic behavior will be significantly changed when the critical non-Hermitian impurity $\lambda g\left(\eta_{x}+i\gamma\eta_{z}\right)$ is applied, that is, any initial state will be forced to evolve to the coalescent state. The corresponding physical observable reflects the property of such a state. This mechanism will be served as the building block to investigate the ODLRO of the nonequilibrium steady state.

In this section, we will demonstrate how does the local impurity can dramatically change the structure of the system. To proceed, we consider $H_{I}$ of Eq. (3) in which the local impurities are applied to a part of lattice sites. An extreme case is that the local dissipation is applied to one lattice site. Compared with the homogeneous case, the $\eta$-pairing states $\{\left|N,l\right\rangle\}$ cannot form an invariant subspace due to the relation of $\left[H_{I},\text{ }H_{0}\right]\neq 0$. Consequently, we cannot judge the EP of the whole system only by the performance of $H_{I}$. However, we can infer the property of $H$ by its counterpart $\overline{H}$ that can be obtained by the transformation $\mathcal{S=}$ $e^{i\theta\eta^{y}}$ with $\tan^{-1}\theta=i\gamma/\lambda$. It represents a rotation in the $\eta^{x}$-$\eta^{z}$ plane around $\eta^{y}$ axis by an complex angle $\theta$ determined by the strength of local pair fluctuation and imaginary field. Notably, the rotation operator is valid at arbitrary $\gamma/\lambda$ unless at the EP ($\gamma/\lambda=\pm 1$) of $H_{I}$, where $\bm{h}_{j}\cdot\bm{\eta}_{j}$ cannot be diagonalized. Applying the transformation $\mathcal{S}$, $H$ is transformed to $\overline{H}=H_{0}+\sqrt{\lambda^{2}-\gamma^{2}}\sum_{j}g_{j}\zeta_{j}^{x}$ where the new set of operators $\zeta_{j}^{\pm,z}=\mathcal{S}^{-1}\eta_{j}^{\pm,z}\mathcal{S}$ also obey the Lie algebra, that is $[\zeta_{j}^{+},$ $\zeta_{j}^{-}]=2\zeta_{j}^{z}$ and $[\zeta_{j}^{z},$ $\zeta_{j}^{\pm}]=\pm\zeta_{j}^{\pm}$. Notice that $\zeta_{j}^{\pm,z}\neq(\zeta_{j}^{\mp,z})^{\dagger}$ owing to the complex angle $\theta$. However, the matrix form of $\overline{H}$ is Hermitian under the biorthogonal basis of $\{\mathcal{S}^{-1}\left|\phi\right\rangle\}$ and $\{\mathcal{S}^{\dagger}\left|\phi\right\rangle\}$ when $\left|\lambda\right|>\left|\gamma\right|$, where $\left\{\left|\phi\right\rangle\right\}$ represents all the eigenstates of $\eta^{z}$. It is worth pointing out that such a transformation solely depends on the ratio of $\gamma/\lambda$ and hence the spectrum is entirely real, even though a nonzero inhomogeneous $g_{j}$ presents. Furthermore, it also indicates that the presence of the local pair fluctuation and imaginary field breaks the SU(2) symmetry of the system but remains the entirely real spectrum without symmetry protection.

Now we investigate the effect of local term on the $\eta$-pairing subspace. The Hermiticity of the matrix form of $\overline{H}$ guarantees the validity of the approximation methods in quantum mechanics. When $\left|\gamma\right|\rightarrow\left|\lambda\right|$, the local term can be treated as a perturbation. With the spirit of the degenerate perturbation theory, the eigenvalues up to the first order are determined by the matrix form of $\sqrt{\lambda^{2}-\gamma^{2}}\sum_{j}g_{j}\zeta_{j}^{x}$ in the subspace spanned by $\{\left|N,l\right\rangle^{\prime}\}$ with $\left|N,l\right\rangle^{\prime}=\mathcal{S}^{-1}\left|N,l\right\rangle$. The corresponding perturbed matrix is referred to as $\overline{H}_{\mathrm{doub}}^{\prime}$, whose elements $\overline{H}_{\mathrm{doub}}^{\prime}$ are given as $\langle\overline{N,l^{\prime}}|\overline{H}_{\mathrm{doub}}^{\prime}\left|N,l\right\rangle^{\prime}=\left(\mathcal{\delta}_{l,l^{\prime}+1}+\delta_{l+1,l^{\prime}}\right)G\sqrt{\lambda^{2}-\gamma^{2}}J_{l}/4N$, where $G=\sum_{j}g_{j}$ and $1/2N$ stems from the translation symmetry of the $\eta$-pairing state. $\{\langle\overline{N,l}|\}$ are biorthogonal left eigenvector in the form of $\{\langle N,l|\mathcal{S}\}$. On the other hand, $H_{I}$ inevitably induce the tunneling from the considered subspace to the other subspace such that the high-order correction to the eigenenergies should be considered. However, the high-order correction term is proportional to the $n$-th power of $\left(\lambda^{2}-\gamma^{2}\right)/U$. If the large $U$ limit and the condition of $\left|\gamma\right|\rightarrow\left|\lambda\right|$ are taken, then one can throw safely high-order correction. Performing the inverse transformation $\mathcal{S}\overline{H}_{\mathrm{doub}}^{\prime}\mathcal{S}^{-1}$, one can obtain the effective Hamiltonian in the doublon invariant subspace

It is a typical non-Hermitian hypercube. When $\left|\lambda\right|=\left|\gamma\right|$, an EP($2N+1$) will be formed no matter which site a local imaginary field is applied to. The corresponding coalescent state is the same with Eq. (13). We stress that such coalescent state is protected by pair interaction during the time evolution.

To verify the conclusion above, we plot the eigenenergies of $H$ as functions of $\gamma/\lambda$ in Fig. 1. Here, we assume that the local imaginary field and pair fluctuation is applied to site $1$ unless stated otherwise. $O_{n}\left(\gamma/\lambda\right)$ is introduced to quantify the similarity between eigenstates $\left|\Phi_{1}\left(\gamma/\lambda\right)\right\rangle$ and $\left|\Phi_{n}\left(\gamma/\lambda\right)\right\rangle$ of $H_{\mathrm{doub}}$, which is defined as

where $\left|\Phi_{1}\right\rangle$ is the eigenstate with lowest eigenenergy. From the Fig. 1, we find that a critical imaginary field cannot only make the degenerate $\eta$-pairing states of $H_{0}$ coalesce but also turn the other degenerate states to a coalescent state and therefore form the multiple high-order EPs of the spectrum. This phenomenon can be understood as follow: for the other eigenstate $\left|\varphi_{0}\right\rangle$ of $H_{0}$, one can also construct an invariant subspace by acting $\eta^{+}$ on such state. The new degenerate subspace belongs to the different eigenvalue of $\eta^{2}$. Taking the same procedures we have done in $H_{\mathrm{doub}}$, a local imaginary field can induce a hypercube-like Hamiltonian as in Eq. (14) which shares the same EP with $H_{\mathrm{doub}}$. The order of EP depends on the degeneracy of eigenenergy of $H_{0}$. This can be demonstrated in Fig. 1. Notice that the system can harbour many EPs with same order that can be identified by the spin symmetry. This result confirms our conclusion on the one hand, and tells us that the cooperation between the local imaginary field and pair fluctuation can accomplish a great task with little effort by clever maneuvers on the other hand. The dramatic change of the eigenstates around EP is the key to understand the following interesting dynamics.

In Yang’s seminar paper Yang (1989), the $\eta$-pairing states are metastable so that they cannot exist stably in the real physical system due to destructive interference from the short-range coherence. In this section, we propose a dynamic scheme to generate a steady state with ODLRO. The scheme is based on the intriguing features of high-order EPs under the influence of a local imaginary field. It can be seen from the previous section that the system spectrum at the EP consists of many types of high-order EPs. As such the system can be decomposed into multiple Jordan blocks owing to the existence of various high-order EPs. In each subspace, an arbitrary initial state will evolve towards the coalescent state and its probability will increase over time in power law according to the order of EP Wang et al. (2016); Yang et al. (2018). Then, a natural question arises: How is the dynamics of such a critical non-Hermitian system? It can be speculated that for an arbitrary initial state, its time evolution depends on the interplay among different types of EP. If the evolution time is long enough, then the highest order EP will determine the final state as the probability in its subspace grows fastest than those in the other subspaces. From this point, the steady state can be generated through dynamical evolution at EP. In the following, we will demonstrate this fascinating behavior through a critical interacting system and investigate the pair correlation of the final evolved state.

Let us consider the critical system with the local imaginary field and pair fluctuation, the Hamiltonian of which can be expressed by setting $j=1$ and $\lambda=\gamma$ of Eq. (3). The scheme is that taking the $\left|N,-N\right\rangle=\left|\mathrm{Vac}\right\rangle$ as an initial state and the final state can be achieved by driving critical non-Hermitian Hamiltonian. According to the aforementioned statement, the effective Hamiltonian $H_{\mathrm{doub}}$ is the key to arrive at the analytical expression of propagator $\mathcal{U=}$ exp$(-iH_{\mathrm{doub}}t)$. For simplicity, we first transform $H_{\mathrm{doub}}$ to a standard Jordan block form (block upper triangular form). Notice that $H_{\mathrm{doub}}$ is a nilpotent matrix with order $2N+1$ and the geometric multiplicity of $\left|\Phi_{\mathrm{c}}\right\rangle$ is $1$. Therefore, $2N$ generalized eigenstates should be introduced to complete this transformation. Such states $\{\left|\Phi_{\mathrm{c}}^{r_{i}}\right\rangle$, $i=1...2N\}$ can be generated from $\left|\Phi_{\mathrm{c}}\right\rangle$ (see Appendix A for more details). Performing transformation $A=\left[\left|\Phi_{\mathrm{c}}\right\rangle\text{, }\left|\Phi_{\mathrm{c}}^{r_{1}}\right\rangle\text{..., }\left|\Phi_{\mathrm{c}}^{r_{2N}}\right\rangle\right]$ on $H_{\mathrm{doub}}$, we can obtain $H_{\mathrm{doub}}^{\mathrm{s}}=A^{-1}H_{\mathrm{doub}}A$ whose matrix element is $\left\langle\overline{N,A_{l}}\right|H_{\mathrm{doub}}^{\mathrm{s}}\left|N,A_{l^{\prime}}\right\rangle=\delta_{l+1,l^{\prime}}g_{1}/2N$ with $\left\{\left|N,A_{l}\right\rangle=A^{-1}\left|N,l\right\rangle\right\}$ and $\left\{\left|\overline{N,A_{l}}\right\rangle=A^{\dagger}\left|N,l\right\rangle\right\}$. Straightforward algebra shows that the propagator in this new frame can be given as

where $h\left(x\right)$ is a Heaviside step function with the form of $h(x)=1$ ($x\geq 0$), and $h(x)=0$ ($x<0$). Notice that the overall phase $e^{-iNU/2}$ is neglected since it does not affect the physical observable. Considering an arbitrary initial state of this subspace, $\sum_{l}c_{l}\left(0\right)\left|N,A_{l}\right\rangle$, the coefficient $c_{m}\left(t\right)$ of the evolved state can be given as

It clearly shows that no matter how the initial state is selected, the coefficient $c_{-N}\left(t\right)$ of evolved state always contains the highest power of time $t$. Consequently, the component $c_{-N}\left(t\right)$ of the evolved state overwhelms the other components ensuring the final state is a coalescent state $\left|N,A_{-N}\right\rangle$ under the Dirac normalization. The different types of the initial state just determines how the total Dirac probability of the evolved state increases over time. Now we return to the original frame, the final state can be given as $\left|\Phi_{\mathrm{c}}\right\rangle=A\left|N,A_{-N}\right\rangle$, which is the coalescent state of $H_{\mathrm{doub}}$ at EP. According to the Appendix B, the correlation of the expected steady state is

It has ODLRO and is independent of the relative distance $r=j-i\neq 0$ between the two operators. This intriguing property is in stark difference to the result obtained in Ref. Tindall et al. (2019), in which the correlation of the steady state decay as the increase of $r$. Furthermore, we investigate the other correlator $\left\langle\Phi_{\mathrm{c}}\right|\eta_{i}^{+}\left|\Phi_{\mathrm{c}}\right\rangle$ that is shown to be zero in the recent experimental and theoretical study Tindall et al. (2019); Kaneko et al. (2019). Straightforward algebra shows that $\left\langle\Phi_{\mathrm{c}}\right|\eta_{i}^{+}\left|\Phi_{\mathrm{c}}\right\rangle=-i/2$. Such quantity is reminiscent of the order parameter $\langle b_{i}^{\dagger}\rangle$ ($b_{i}^{\dagger}$ represents the boson creation operator) that describes the quantum phase transition from the Mott insulator to superfluid in the Bose-Hubbard model. The nonzero value may imply some important physical property and give further insight into the non-Hermitian Hubbard model.

To check this understanding, we perform a numerical simulation and present the results in Fig. 2 with an initial state being $\left|\Psi\left(0\right)\right\rangle=\left|N,-N\right\rangle$. We inspect the time dependence of two correlators $|\left\langle\Psi\left(t\right)\right|\eta_{i}^{+}\left|\Psi\left(t\right)\right\rangle|$ and $\left\langle\Psi\left(t\right)\right|\eta_{i}^{+}\eta_{j}^{-}\left|\Psi\left(t\right)\right\rangle$. It can be shown that the final values of the two correlators are $1/2$, and $1/4$, which agree with our analytical results. The number of local imaginary fields only changes the relaxation time $t_{r}$ and does not change the final value of the correlator. The systems that realize the proposed dynamic scheme are all with the reach of ongoing experiments. Hence, such a scheme offers a unique perspective to generate a steady superconducting-like state in a variety of materials.

In summary, we have investigated some general aspects of the non-Hermitian Hubbard model. We have shown that the presence of the local imaginary field and pair fluctuation can drastically change its microscopic behavior and manifest a variety of collective and cooperative phenomena at the macroscopic level. Specifically, the whole real spectrum of such a non-Hermitian interacting system can be achieved in a wide range of parameters even though a local imaginary field is applied. At EP, the interplay between the local imaginary field and on-site pair fluctuation can lead to a coalescent state with the geometric multiplicity of $1$, which is protected by the on-site pair interaction. $\eta$-pairing symmetry plays the key role to understand the formation of the high-order EP. Comparing with the other schemes in Refs. Tindall et al. (2019); Kaneko et al. (2019) that produce a superconducting-like state, the critical strongly correlated system in our scheme favors superconductivity on a long time scale. The corresponding steady state with ODLRO can be generated through time evolution from an arbitrary initial state. The realization of this scheme does not depend on the location of the local imaginary filed and pair fluctuation. Hence, this scheme does open up the possibility of generating a steady superconducting-like state with a long coherence time in a variety of experimental platforms.