Nonequilibrium phases and phase transitions of the XY-model

We consider an $XY$-spin chain of $N$ sites (labeled by $r$), exchange coupling $J$, and coupled to a transverse field $h$. At its ends, i.e., at $r=1$ and $r=N$, the chain is coupled to magnetic reservoirs which are kept at zero temperature ($T=0$). The Hamiltonian of the chain is given by

where $\sigma_{r}^{x,y,z}$ are the Pauli matrices at site $r$, and $\gamma$ controls the anisotropy. The total Hamiltonian is given by

where ${\cal H}_{l}$ and ${\cal H}_{\text{C-}l}$, with $l=\text{L,R}$, are respectively the Hamiltonians of the reservoirs and the system-reservoir coupling terms. In the following, we assume that the reservoirs possess bandwidths which are entirely determined by magnetic potential $\mu_{l}$ ($l=\text{L,R}$) and which are much larger than the energy scales that characterize the chain.

In the wide-band limit, results become independent of the details of ${\cal H}_{l}$ and ${\cal H}_{\text{C-}l}$. For concreteness, we take the reservoirs to be isotropic $XY-$chains, i.e.

with $l=\text{L},\text{R}$, and we have defined $\Omega_{\text{L}}\equiv\left\{-\infty,\ldots,0\right\}$, $\Omega_{\text{R}}\equiv\left\{N+1,\ldots,\infty\right\}$. Initially, the reservoirs are in an equilibrium Gibbs state, $\rho_{l}=e^{-\beta(\mathcal{H}_{l}-\mu_{l}M_{l})}$, where $M_{l}=\sum_{r_{l}\in\Omega_{l}}\sigma_{r_{l}}^{z}$ is the reservoir magnetization (which is a good quantum number in the absence of system-reservoir coupling, i.e., $\left[{\cal H}_{l},M_{l}\right]=0$). The average value of $M_{l}$ is set by the magnetic potential $\mu_{l}$. For finite $\mu_{l}$ these are non-Markovian reservoirs, with power-law decaying correlations, and a set of gapless magnetic excitations within an energy bandwidth $J_{l}\gg J,h$. The chain-reservoir coupling Hamiltonians are

with $\left(r_{\text{L}}\right)_{\text{C}}=1$, $\left(r\right)_{\text{L}}=0$, $\left(r_{\text{R}}\right)_{\text{C}}=N$, and $\left(r\right)_{\text{R}}=N+1$. A sketch of this system is shown in Fig.1(a).

The full Hamiltonian, ${\cal H}$, can be represented in terms of fermions via the so-called Jordan-Wigner (JW) mapping(Lieb et al., 1961a), $\sigma_{r}^{+}=e^{i\pi\sum_{r^{\prime}=1}^{r-1}\hat{c}_{r^{\prime}}^{\dagger}\hat{c}_{r^{\prime}}}c_{r}^{\dagger}$, where $\hat{c}_{r}^{\dagger}/\hat{c}_{r}$ creates/annihilates a spinless fermion at site $r$. The JW-transformed system corresponds to a Kitaev chain(Kitaev, 2001) in contact with two metallic reservoirs of spinless fermions at chemical potentials $\mu_{\text{L},\text{R}}$, i.e.

where $J\gamma$ defines the superconducting coupling strength and $h$ plays the role a potential applied on the chain. A sketch of this system is shown in Fig.1(b). In equilibrium, topologically non-trivial phases of the Kitaev chain correspond to magnetically ordered phases of the original XY-spin model, whereas the topologically trivial cases correspond to disordered phases. With a magnetic bias, the transfer of spin excitations between the reservoirs was studied rather extensively (see, e.g., Refs. Meier and Loss, 2003; Nakata et al., 2017), also considering transport signatures of the topological phase in short junctions.Hoffman et al. (2018); Shen et al. (2020) Here, however, we will be mostly concerned with the bulk properties at $N\to\infty$, and after the reservoirs have been traced out.

As the JW-transformed Hamiltonian is quadratic in its fermionic degrees of freedom, the non-equilibrium system admits an exact solution in terms of single-particle quantities. In the following, we employ the non-equilibrium Green’s function formalism to compute correlation functions and related observables. The procedure is described in the Supplemental Material of Ref. [Puel et al., 2019] and is briefly summarized here for convenience.

We start by defining the Nambu vector, $\hat{\boldsymbol{\Psi}}^{\dagger}=\left(\hat{c}_{1}^{\dagger},\ldots,\hat{c}_{N}^{\dagger},\hat{c}_{1},\ldots,\hat{c}_{N}\right)$, and the retarded, advanced, and Keldysh components of the Green’s function, given by

Using this notation, the Hamiltonians for the right, left reservoirs and for the chain are given by ${\cal H}_{l}=\frac{1}{2}\hat{\boldsymbol{\Psi}}^{\dagger}\boldsymbol{H}_{l}\hat{\boldsymbol{\Psi}}$, with $l=\text{L,R,C}$. For the chain, $\boldsymbol{H}_{\text{C}}$ is a $2N\times 2N$ Hermitian matrix respecting particle-hole symmetry, i.e., $\boldsymbol{S}^{-1}\boldsymbol{H}_{\text{C}}^{T}\boldsymbol{S}=-\boldsymbol{H}_{\text{C}}$ where $\boldsymbol{S}=\tau^{x}\otimes\boldsymbol{1}_{N\times N}$ and $\tau^{x}$ interchanges particle and hole spaces. Similar definitions apply to the degrees of freedom of the right and left reservoirs. The bare retarded and advanced Green’s functions, in the absence of chain-reservoir couplings, are simply given by $\boldsymbol{G}_{l,0}^{R/A}\left(\omega\right)=\left(\omega-\boldsymbol{H}_{l}\pm i\eta\right)^{-1}$. Re-writing the chain-reservoir coupling in the same notation, ${\cal H}_{\text{C-}l}=\frac{1}{2}\left(\hat{\boldsymbol{\Psi}}^{\dagger}_{(l)}\boldsymbol{T}\hat{\boldsymbol{\Psi}}+\hat{\boldsymbol{\Psi}}^{\dagger}\boldsymbol{T}^{\dagger}\hat{\boldsymbol{\Psi}}_{(l)}\right)$. The self-energy of the chain, induced by tracing out the reservoirs, are $\boldsymbol{\Sigma}^{R/A/K}=\sum_{l=\text{L,R}}\boldsymbol{\Sigma}^{R/A/K}_{l}$ where

are the contributions of reservoir $l$, which obey equilibrium fluctuation-dissipation relations, and where $n_{\text{F},l}(\omega)=(e^{\beta_{l}(\omega-\mu_{l})}+1)^{-1}$ is the Fermi-function with chemical potential $\mu_{l}$ and inverse temperature $\beta_{l}$. The chain steady-state Green’s functions are obtained from the Dyson’s equation,

As mentioned above, we consider the case where the bandwidths of the reservoirs, $J_{l=\text{L},\text{R}}$ are much larger than the other energy scales. In this wide band limit, the coupling to reservoir $l$ is completely determined by the hybridization energy scale $\Gamma_{l}=\pi J_{\text{C-}l}^{2}D_{l}$. Here, $D_{l}$ is the reservoir’s constant-local density of states. In practice, the wide band limit yields a frequency independent retarded self-energy, $\boldsymbol{\Sigma}^{R}_{l}=i(\boldsymbol{\gamma}_{l}+\bar{\boldsymbol{\gamma}}_{l})$, which substantially simplifies subsequent calculations, with $\boldsymbol{\gamma}_{l}=\Gamma_{l}\left|r_{l}\right\rangle\left\langle r_{l}\right|$ and $\bar{\boldsymbol{\gamma}}_{l}=\Gamma_{l}\left|\bar{r}_{l}\right\rangle\left\langle\bar{r}_{l}\right|$, and where $\left|r\right\rangle$ and $\left|\bar{r}\right\rangle\equiv\boldsymbol{S}\left|r\right\rangle$ are single-particle and hole states.

In this case, it is convenient to define the non-Hermitian single-particle operator

which we assume to be diagonalizable, possessing right and left eigenvectors $\left|\alpha\right\rangle$ and $\left\langle\tilde{\alpha}\right|$, and associated eigenvalues $\lambda_{\alpha}$. In terms of these quantities, the retarded Green’s function is given by

and the Keldysh Green’s function becomes

Steady-state observables can be obtained from the single-particle correlation function matrix, $\boldsymbol{\chi}\equiv\langle\hat{\boldsymbol{\Psi}}\hat{\boldsymbol{\Psi}}^{\dagger}\rangle$, which is obtained from the Keldysh Green’s function

The explicit form of $\boldsymbol{\chi}$ after performing the integration over frequencies is provided in Eq.(41). As the model is quadratic, $\boldsymbol{\chi}$ encodes all the information about the reduced density matrix of the chain, $\hat{\rho}_{\text{C}}=\text{tr}_{\text{L,R}}[\hat{\rho}]$. This quantity can itself be expressed as the exponential of a quadratic operator, i.e. $\hat{\rho}_{\text{C}}=\text{e}^{\hat{\Omega}_{\text{C}}}/Z$, where $Z=\text{tr}\left[\text{e}^{\hat{\Omega}_{\text{C}}}\right]$ and $\hat{\Omega}_{\text{C}}=\frac{1}{2}\hat{\boldsymbol{\Psi}}^{\dagger}\boldsymbol{\Omega}_{\text{C}}\hat{\boldsymbol{\Psi}}$ with $\boldsymbol{\Omega}_{\text{C}}$ being a $2N\times 2N$ matrix respecting the particle-hole symmetry conditions. $\hat{\Omega}_{\text{C}}$ is related to the single-particle density matrix via

This relation allows the calculation of mean values of quadratic observables, $\hat{O}=\frac{1}{2}\hat{\boldsymbol{\Psi}}^{\dagger}\boldsymbol{O}\hat{\boldsymbol{\Psi}}$, defined by the Hermitian and particle-hole symmetric matrix $\boldsymbol{O}$,

as well as all higher-order correlation functions.

The system in equilibrium is more conveniently studied without considering the couplings to the leads and assuming periodic boundary conditions. After performing the JW transformation, the Hamiltonian of the translation-invariant chain is diagonalized in the momentum representation by a suitable Bogoliubov transformation, i.e., ${\cal H}_{C}=\sum_{k}\varepsilon_{k}(\hat{\gamma}_{k}^{\dagger}\hat{\gamma}_{k}-1/2)$, where the operators $(\hat{\gamma}_{k},\hat{\gamma}_{-k}^{\dagger})^{T}=e^{i\theta_{k}\sigma_{x}}(\hat{c}_{k},\hat{c}_{-k}^{\dagger})^{T}$ describe excitations of energy

and $\sin\left(2\theta_{k}\right)=-2J\gamma\sin(k)/\varepsilon_{k}$.

The ground state is characterized by a vanishing number of Bogoliubov excitation, i.e., $n_{k}=0$ where

For $\left|h/J\right|<1$, the ground-state is topologically non-trivial with positive and negative anisotropies ($\gamma>0$ or $<0$) corresponding to opposite signs of the topological invariant, separated by a critical gapless state for $\gamma=0$. At $\left|h/J\right|=1$, the spectral gap vanishes and the system transitions into a topologically trivial phase at large $h$.

A similar phase diagram is obtained in terms of the original spin degrees of freedom. Fig. 2(a) illustrates the zero-temperature phase diagram of the equilibrium $XY$-model. For $\gamma>0$, the systems is magnetically ordered along the $x$ direction under a weak transverse field, $h$, and possesses a finite magnetization(Barouch and McCoy, 1971) $\phi\equiv\lim_{h_{x}\to 0}\lim_{L\to\infty}\frac{1}{N}\sum_{r}\left\langle\sigma^{x}_{r}\right\rangle\neq 0$, where $h_{x}$ is a symmetry-breaking magnetic field along the $x$ direction. A negative anisotropy, i.e., $\gamma<0$, yields a non-vanishing magnetization along the $y$ direction, whereas at $\gamma=0$ the system is critical and isotropic for $\left|h/J\right|<1$. As the ordered phases for $\gamma>0$ and $<0$ are equivalent to each other and related via a simple rotation, only $\gamma>0$ is considered in the subsequent analysis. It is worth recalling that the special cases $\gamma=\pm 1$ correspond to the transverse field Ising model. A strong $h$ drives the magnetic phase through a second order phase transition into a phase of vanishing magnetization, i.e., a phase with $\phi=0$. Near the transition, for $\left|h/J\right|<1$, the magnetization behaves as $\phi\simeq\sqrt{\frac{2}{1+|\gamma|}}\left\{\gamma^{2}\left[1-\left(h/J\right)^{2}\right]\right\}^{1/8}$ Barouch and McCoy (1971).

The computation of the order parameter, $\phi$, directly from the above definition is not possible via the JW mapping. Instead one considers the two-points correlation functions ($\alpha=x,y,z$):

For disordered phases in equilibrium, these correlators are expected to show either exponential (EXP) or power-law (PL) decay depending on whether the system is gapped or gapless. In the ordered phase the system has long-range order (LRO) correlations, e.g., for $\gamma>0$,

where $\xi$ is the characteristic correlation length and $A$ a numeric coefficient. This expression allows one to obtain $\phi$ from the correlation function $\mathbb{C}_{r,r^{\prime}}^{xx}$, which in turn can be computed in terms of a Toeplitz determinant (Lieb et al., 1961b; Sachdev, 1996). In equilibrium, all correlation functions, except $\mathbb{C}_{r,r^{\prime}}^{xx}$, either vanish or decay exponentially with $\left|r-r^{\prime}\right|$.

For an open system connected to demagnetized baths, i.e., $\mu_{L,R}=0$, the same equilibrium bulk properties as those for the closed system are found. It is natural to expect that bulk properties pertain for distances greater than $\xi$ away from the leads. Our calculation of fermionic observables follows that described in Sec.II. The calculation of the spin-spin correlation functions are similar to those for the translation-invariant system and is given in Appendix A in terms of the single-particle correlation matrix $\boldsymbol{\chi}$.

The non-equilibrium features of the XY-model with Markovian reservoirs were first reported in Refs. [Prosen and Pižorn, 2008; Banchi et al., 2014]. This limit can be recovered from the present model by taking $\left|\mu_{\text{R}}\right|$ or $\left|\mu_{\text{L}}\right|\rightarrow\infty$(Ribeiro and Vieira, 2015a). The steady-state phase diagram in that limit possesses two distinct phases characterized respectively by an exponential decay of all correlation functions with distance and by an algebraic decay of $\mathbb{C}^{zz}_{r,r^{\prime}}$ concomitant with a strong sensitivity to small variations of some control parameters(Prosen and Pižorn, 2008). Fig. 2(b) depicts the phase diagram in this Markovian limit and with its sensitive (white) and normal (light-yellow) phases. The dark yellow lines mark critical phase boundaries.

The quasiparticle dispersion relation, given by Eq.(19), is shown in 3(a) and (c). The algebraic correlations are associated to the presence of an inflection point in the quasiparticle dispersion which appears for $h/J\leq|1-\gamma^{2}|$. In the normal region, the extrema of the energy are $m_{1}=\varepsilon_{k=0}=2J\left|1+h/J\right|$ and $m_{2}=\varepsilon_{k=\pi}=2J\left|1-h/J\right|$, while for the sensitive region $m_{3}=2J\left|\gamma\right|\sqrt{1+\left(h/J\right)^{2}\left(\gamma^{2}-1\right)^{-1}}$ becomes a global extremum. Figs. 3(b) and 3(d) show the spectrum of the non-Hermitian single-particle operator $\boldsymbol{K}$ (see Eq. (13)) for both normal and sensitive phases. It turns out that the imaginary part of the eigenvalues scales with the inverse system size, $\text{Im}\lambda_{\alpha}\propto N^{-1}$. (Medvedyeva and Kehrein, 2014) This is a reflection of the fact that for the chain degrees of freedom, the dissipative effects of the boundary become less important with increasing system size. A key feature of the sensitive region is that, for energies (i.e. $\text{Re}(\lambda_{\alpha})$) where four momenta can propagate, the spectrum does not converge to a line with increasing system size, but becomes scattered within a finite area (Prosen and Pižorn, 2008). These effects are independent of the Markovian nature of the reservoirs and remain for the non-Markovian case as the operator $\boldsymbol{K}$ does not depend on the chemical potential of the leads. Thus, as explicitly shown below, the normal-sensitive transition, reported in Refs. [Prosen and Pižorn, 2008; Banchi et al., 2014] for the Markovian case, also occurs for finite values of $\mu_{\text{L}}$ and $\mu_{\text{R}}$.

We now present the phase diagram of non-equilibrium $XY$ model and show that the energy current passing through the chain can be used to discriminate between the different phases.

Conservation of energy implies that the steady-state energy current is equal across any cross section along the chain and can be obtained from $\boldsymbol{\chi}$ as

where $\boldsymbol{J}_{r}$ is the single-particle current operator at link $(r,r+1)$ which is explicitly given in Appendix A.

We have previously discussed the steady-state energy current in a non-Markovian setting for the particular case of the transverse field Ising model, i.e., $\left|\gamma\right|=1$, in Ref. [Puel et al., 2019]. The non-equilibrium phase diagram, as a function of $\mu_{\text{L}}$ and $\mu_{\text{R}}$, in the normal phase, is qualitatively similar to the Ising case and is reproduced in Fig. 4(a). Two of the phases which arise near $\mu_{\text{L}}=\mu_{\text{R}}$, do not support energy transport,i.e., $\mathcal{J}_{e}=0$: the ordered phase (O) and the non-conducting phases (NC). Other phases may be further characterized in terms of their energy conductance, i.e., ${\cal G}_{\text{L}}\equiv\partial_{\mu_{\text{L}}}\mathcal{J}_{e}$ and ${\cal G}_{\text{R}}\equiv\partial_{\mu_{\text{R}}}\mathcal{J}_{e}$. We refer to current-saturated (CS) the phases with $\mathcal{J}_{e}\neq 0$ and ${\cal G}_{\text{R}}={\cal G}_{\text{R}}=0$. They arise when one of the reservoirs chemical potentials is larger than $m_{1}$ while the other lies inside the quasiparticle excitation gap. The conducting phase (C) is characterized by a non-zero conductance, i.e., ${\cal G}_{\text{L}}\neq 0$ and/or ${\cal G}_{\text{R}}\neq 0$, arising whenever at least one of the chemical potentials lies within the quasiparticle excitations band, i.e., $\left|\mu_{\text{L}}\right|$ and/or $\left|\mu_{\text{R}}\right|\in(m_{1},m_{2})$.

Fig. 4(b) depicts the phase diagram for a generic $XY$ chain. Besides the phases found for $\gamma=1$, an analysis of the occupation numbers (see next subsection) shows that some regions acquire a noise-like behavior. These phases, similar to the sensitive regions of the Markovian case, are labelled NC^∗, CS^∗, and C^∗.

In Figs. 4(c) and 4(d) we show the current $\mathcal{J}_{e}$ and conductance ${\cal G}_{\text{L}}$ for a fixed $\mu_{\text{R}}$ represented by the red-dashed lines in the phase diagrams. In the sensitive region, the C phase is crossed by the transition line at $\mu_{L}=m_{2}$, where the conductance becomes non-analytic. It is worth noting that the noise-like behavior found in the occupation numbers does not appear in the current of energy.

In terms of the JW fermions, the present analysis is similar to that of a transport across a tight binding model in the sense that when the chemical potentials cross the dispersion relation non-analytic properties of the current appears. However, the existence of anomalous terms in the fermionic Hamiltonian inhibits a closer comparison with charge transport.

In the current carrying steady-state regime of a Fermi gas, fluctuations in the number of particles were shown to be intimately related to the entanglement entropy of a subsystem (Klich and Levitov, 2009; Song et al., 2010, 2011, 2012). In analogy, the occupation number of the Bogoliubov excitations, given by Eq.(20), can be used to describe the properties of the asymptotic steady-state away from the boundaries. In the open system setting, $n_{k}$ can be approximated by numerically computing the Fourier transform $\boldsymbol{\chi}_{k}=\sum_{r\in\Omega}e^{-ik(r-r_{0})}\boldsymbol{\chi}_{r,r_{0}}$, where $\Omega=\left\{r:L/4<r<3L/4\right\}$ and $r_{0}=L/2$, followed by a Bogoliubov transformation. In equilibrium, $\mu_{L}=\mu_{R}=0$, $n_{k}\simeq 0$ as expected, while in a generic out-of-equilibrium situation $n_{k}\neq 0$.

For the Ising model, it was shown that in the CS and NC phases the $n_{k}$ is a continuous function of $k$, while in the C phase it has discontinuities depending on the reservoir’s chemical potentials (Puel et al., 2019). These discontinuities happen at the momenta $\pm k_{l=\text{L,R}}$ where the chemical potential $\mu_{l=\text{L,R}}$ cross the dispersion relation, see Fig. 5(a). These results extent straightforwardly to the XY-model in the normal region, see Fig. 5(c)-(f). Within the O phase the system behaves as in equilibrium, i.e. $n_{k}\simeq 0$.

In the sensitive region there may be two absolute values of momenta, labelled $\pm k_{l=\text{L,R}}$ and $\pm k^{\prime}_{l=\text{L,R}}$, for which each chemical potential crosses the dispersion relation, as illustrated in Fig. 5(b). Interestingly, we find that $n_{k}$ has an intrinsic noise in the sensitive region, see Figs. 5(g)-(j). The noise appears in phases NC^∗, C^∗, and CS^∗, for $\left|k\right|>k_{m_{2}}$, where $\varepsilon_{k_{m_{2}}}=m_{2}$ (see Fig. 3(c)). In Appendix B we check that the magnitude of the noise in $n_{k}$ does not diminishes with increasing system sizes. Curiously, the noise vanishes along the line $\mu_{\text{L}}=-\mu_{\text{R}}$, as well as within phases C and CS crossed by this line, as shown in Figs. 5(k) and (l), and studied in detail in Appendix B.

Note that $n_{k}$ is asymmetric upon changing $k\to-k$ for all conducting phases as required to maintain a net energy flow through the chain, as $\varepsilon(k)=\varepsilon(-k)$. Figs. 5(c)-5(l) illustrates this feature by showing a larger value of the hybridization, that yield a larger current of energy and consequently to a more asymmetric $n_{k}$.

As aforementioned, for $\gamma>1$, the magnetization along the $x$ direction, $\phi$, is a good order parameter for the broken-symmetry equilibrium phase. In the open system $\phi$ can still be used as order parameter. However, by changing the chemical potential of the, say, left reservoir, $\phi$ drops to zero discontinuously as soon as the system reaches the disordered phase. Interestingly, this transition shows a mixed-order behavior where the discontinuity of $\phi$ is accompanied by a divergence of the correlation length (Puel et al., 2019). The divergence occurs as $\xi\propto\left|\mu_{\text{L}}\pm m_{i}\right|^{-\nu}$ when approaching the critical point from the disordered phase ($i=2,3$ depending on the values of $\gamma$ and $h$). The critical exponent is $\nu=1/2$, except for special values of the parameters (see Appendix C).

We have numerically verified that this behavior survives away from $\gamma=1$, with the same type of dependence of the order parameter and correlation length. In particular, the mixed-order transition with $\nu=1/2$ is also present in the sensitive region, e.g., at the transitions between O and C^∗ phases in the phase diagram of Fig. 4(b). We show in panels (c) and (d) of Fig. 6 an example of the numerical analysis of order parameter and correlation length in the sensitive region. The main difference compared to the normal region is that here finite size effects are much stronger, which is consistent with the sensitive dependence on $N$ of the the rapidity spectrum and occupation numbers, see Figs. 3 and 5.

To derive the value of the critical exponent $\nu$, we consider the explicit form of the correlation function in terms of a Toepliz determinant:(Lieb et al., 1961b; Sachdev, 1996)

where $D_{n}$ is given by:

The asymptotic dependence can be obtained from Szego’s lemma, leading to the following expression for the correlation length:

Here, the difference from the standard treatment of the transverse-field Ising chain(Lieb et al., 1961b; Sachdev, 1996) is simply that the occupation numbers $n_{k}$ are kept generic, thus are allowed to assume any non-equilibrium distribution induced by the external reservoirs. For example, Eq. (26) takes into account that in general $n_{k}\neq n_{-k}$, as shown by Fig. 5 with a large hybridization energy. By substituting the Fermi distribution, Eqs. (25) and (26) recover the known equilibrium expressions at finite temperature.Sachdev (1996, 2011)

Applying Eq. (26) to the critical point, we first assume as in Fig. 5(a) that the minimum of the quasiparticle dispersion occurs at $k=\pi$. Furthermore, if $\mu_{\rm R}$ is inside the gap the critical point is at $\mu_{\rm L}=m_{2}$. Close to the critical point, the only occupied states are in a small range $k\in[\pi-\Delta k_{\rm L},\pi+\Delta k_{\rm L}]$ around the minimum of $\varepsilon_{k}$, thus we can approximate the quasiparticle dispersion as parabolic giving $\Delta k_{L}\propto\sqrt{\mu_{\rm L}-m_{2}}$. This dependence of $\Delta k_{L}$ is directly related to the critical exponent $\nu=1/2$. More precisely, $\Delta k_{L}$ close to the critical point is given by:

and we can set $n_{k}\simeq n_{k=\pi}$ in the small integration interval of Eq. (26), leading to:

This expression clearly shows how the divergence of $\xi$ is due to the shrinking of the region of nonzero occupation. We show in Fig. 8 that this theory is accurate, by a direct comparison to the numerical results.

After having clarified the origin of the critical exponent $\nu=1/2$, it is interesting to compare the behavior of the open chain to the temperature dependence of the equilibrium system. In the latter case, an ordered phase is only allowed at zero temperature and the order disappears at any arbitrarily small temperature $T>0$. The sudden disappearance of the ordered state is related to the presence of thermal excitations, and is analogous to the vanishing of $\phi$ induced by the non-equilibrium chemical potentials, as soon as either $\mu_{\rm L}$ or $\mu_{\rm R}$ overcomes the gap. Furthermore, similarly to the non-equilibrium system, the correlation length diverges when $T\to 0$. As it turns out, in the low-temperature limit, $\xi$ can be related in a simple way to the density of excitations $\rho=\int_{0}^{2\pi}\frac{dk}{2\pi}n_{k}$:

This expression follows immediately from Eq. (26), since $n_{k}=n_{-k}\ll 1$ when $T\to 0$, and can also be understood by a simple argument in terms of a dilute gas of domain-wall excitations.Sachdev (2011)

A naive application of Eq. (29) to the non-equilibrium system is shown in Fig. 8. Although Eq. (29) predicts the correct critical exponent $\nu=1/2$, there is a clear disagreement with the numerical results. This failure of Eq. (29) can be explained from the non-vanishing value of $n_{k}$ around the minimum of $\varepsilon_{k}$ (say, $k=\pi$): in both the low-temperature limit and the non-equilibrium system we have $\rho\to 0$ at the critical point, which results in a diverging correlation length. However, in the first case we have $n_{k}\to 0$ while for the non-equilibrium system $n_{\pi}$ remains finite, and the vanishing of $\rho$ is due to the shrinking of $\Delta k_{L}$. Instead of the density $\rho$, we can consider a ‘modified’ density:

which follows naturally from Eq. (26) by requiring that a relation similar to the equilibrium system at low temperature is satisfied, $\xi^{-1}=2\rho_{\xi}$. It is easy to see that close the critical point we have $\rho_{\xi}\simeq-(\log\left[1-2n_{\pi}\right]/2n_{\pi})\rho$, which differs from $\rho$ by a nontrivial multiplicative factor. An interesting exception, discussed more extensively in Appendix C, occurs for $J=h/2$, when $n_{\pi}=0$ and the relation between $\xi$ and $\rho$ is Eq. (29), as in equilibrium. At $J=h/2$ the vanishing of $n_{k}$ also affects the value of the critical exponent, which is $\nu=5/2$ instead of 1/2.

The above arguments can be adapted to other parameter regimes. In particular, the discussion is almost unchanged for $h/J<\min[0,\gamma^{2}-1]$, when the dispersion minimum is at $k=0$. Instead, the treatment of the sensitive phase with $|\gamma|<1$ is more delicate. Firstly, as shown in Fig. 5(b), the dispersion is characterized by two minima instead of one. More importantly, $n_{k}$ appears to be highly pathological when $N\to\infty$, when the occupation numbers undergo wild oscillations. Despite these differences, numerical evaluation of the critical exponent still gives $\nu=1/2$, thus we conjecture that a suitable average of $n_{k}$ is well-defined:

Then, the critical exponent would be determined through Eq. (28) in a way analogous to the regular case, i.e., the critical exponent would correspond to the shrinking of the occupied regions around the minima, explaining the persistence of $\nu=1/2$ in this phase.

As summarized in Table 1, $\mathbb{C}_{r,r^{\prime}}^{zz}$ displays a power-law decay in several of the allowed phases. We focus here on the the non-sensitive region, where this behavior can be understood from the well-known power-law decay of density-density correlations of non-interacting fermions. In fact, through the fermionic mapping, $\mathbb{C}_{r,r^{\prime}}^{zz}$ is equivalent to a density-density correlation function:

In the C phase, at least one of the reservoirs has its chemical potential within the range of quasiparticle energy spectrum. Then, the correlation function is expected to have a power-lay decay with oscillating character, similar to a simple 1D Fermi gas where it decays as $|r-r^{\prime}|^{-2}$ and oscillates with wavevector $2k_{F}$, with $k_{F}$ the Fermi wavevector (see, e.g., Ref. Castin, 2007).

In our case, we express the correlation function through the Bogolubov excitations $\hat{\gamma}_{k}$ of the transnational invariant system (which is appropriate in the thermodynamic limit). The corresponding occupation numbers are defined in Eq. (20) and give:

If for simplicity we assume $n_{k}\simeq n_{-k}$ (which is justified in the limit of vanishing hybridization energy $\Gamma$) the occupation numbers have discontinuities at $k=\pm k_{i}$, induced by the left and right reservoirs ($i=L,R$). When $k_{i}r\gg 1$, we can extract the leading contribution to Eq. (IV.2)induced by the discontinuous jumps $\Delta n_{k_{i}}$ of $n_{k}$, defined by $\partial_{k}n_{k}=\sum_{i=L,R}\Delta n_{k_{i}}\left[\delta(k+k_{i})-\delta(k-k_{i})\right]$:

which simplifies to:

when there is a single Fermi surface (here, induced by $i=R$). The above expressions display the expected $1/r^{2}$ decay and oscillatory dependence. As shown in Fig. 9, we find a good agreement between Eq. (35) and the numerical results.

To better characterize the oscillatory dependence, we have also studied the Fourier transform:

which is shown in Fig. 10 for two representative cases. With a single Fermi surface (left panels) we find dominant non-analytic features at $k=0,2k_{R}$, in agreement with Eq. (35). We also find smaller discontinuities in $\partial_{k}\mathbb{C}_{k}^{zz}$ at higher harmonics, $k=4k_{R},6k_{R}$, which are not captured by the leading-order approximation Eq. (35). With two Fermi surfaces (right panels) we find the expected singularities at $k=0,2k_{L,R}$. However, there are additional features at $\partial_{k}\mathbb{C}_{k}^{zz}$ at $k=k_{L}\pm k_{R}$, which are in agreement with Eq. (IV.2). As seen there, the correlation function is not simply a sum of $i=L,R$ contributions, but involves interference terms between the two Fermi surfaces.

Finally, we comment on the power-law dependence of $\mathbb{C}_{r,r^{\prime}}^{zz}$ in the sensitive region. Away from the ordered phase we find a power-law decay $|r-r^{\prime}|^{-s}$ where, however, the exponent is generally different from $s=2$ (we often find $s<2$) and depends on system parameters. This behavior is most likely related to the singular nature of $n_{k}$, which from our numerical evidence is characterized by a complex pattern of closely spaced discontinuities (see, e.g., Fig. 5). Such discontinuities will contribute to the square parenthesis of Eq. (IV.2) in a way difficult to compute explicitly (the summation index $i$ should become a continuous parameter) and might be able to modify the exponent $s$. This interpretation is confirmed by the survival of the power-lay decay in the NC^∗ and CS^∗ regions, where the chemical potentials $\mu_{\rm L,R}$ do not cross the quasiparticle bands, thus an exponential decay might be expected. Instead, Fig. 5(g) and (i) show that a discontinuous dependence of $n_{k}$ can be found in these regions as well, in agreement with the observed power-law dependence of $\mathbb{C}_{r,r^{\prime}}^{zz}$. Finally, the simple discontinuities of Fig. 5(l) result in the regular value $s=2$.