Non-Hermitian quantum gases: a platform for imaginary time crystals

One of the most important applications of quantum mechanics is the thermodynamic description of quantum gases. Despite the fundamental importance of this topic, a comprehensive description of the thermodynamic properties of non-Hermitian quantum gases is still lacking. Here, we investigate the properties of bosonic and fermionic non-Hermitian systems at finite temperatures. We show that non-Hermitian systems exihibit oscillations both in temperature and imaginary time. As such, they can be a possible platform to realize an imaginary time crystal (iTC) phase. The Hatano-Nelson model is identified as a simple lattice model to reveal this effect. In addition, we show that the conditions for the iTC to be manifest are the same as the conditions for the presence of disorder points, where the correlation functions show oscillating behavior. This analysis makes clear that our realization of iTC is effectively a way to filter one specific Matsubara mode. In this realization, the Matsubara frequency, that enters as a mathematical tool to compute correlation functions for finite temperatures, can be measured experimentally.


Introduction
The evolution of a system in imaginary time τ is a long known prescription in quantum field theory to compute the partition function [1]. In this formalism, periodic boundary conditions lead to discrete sets of frequencies, the Matsubara frequencies ω n . Despite their pivotal importance, these frequencies are believed to be just tools to compute correlation functions. Observing the structure of the Green's function G, one concludes that there will be no poles associated with ω n , but only to the real frequencies, with a small imaginary part related to the causal structure of G. A natural question emerges then when one considers systems with sizable values of the imaginary part of theirs energies, as it occurs for non-Hermitian systems.
The simplest kind of non-Hermitian system is a non-Hermitian quantum gas. A quantum gas is a non-interacting system, such that its HamiltonianĤ can be decomposed as the sum of Hamiltonians of some quantum numbers m, H = ⊕ mĥm . The thermodynamic properties of a non-Hermitian ideal quantum gas were not explored so far; yet, in the presence of a pseudo-Hermitian symmetry, a biorthogonal thermodynamic description is available in the literature [18]. Although the thermodynamic limit is not really achievable for large OBC systems due to the NHSE, a finite-size thermodynamic analysis is still consistent and one can obtain results in the thermodynamic limit for PBC and for the surrogate Hamiltonian (SH). The latter consists of an analytical continuation of the Bloch Hamiltonian that reproduces the bulk spectrum of the system with OBC [7,9,19].
Here, we investigate the thermodynamic behavior of non-Hermitian quantum gases at finite temperature. We find that for a range of intermediary temperatures, this system exhibits oscillations in both β = 1/(k B T ) and imaginary time τ . This is precisely a footprint of the imaginary time crystal (iTC) phase conjectured by Wilczek in his original paper on time crystals [20]. The conditions for the existence of this phase are the same as the ones 2 for the presence of disorder points [21][22][23][24][25][26][27]-critical phases at which the correlation function has a modulation, together with the exponential decay. We apply our results to the Hatano-Nelson model [3,28] and show how one can observe these oscillating phases in both space and imaginary time.

Thermodynamics of the Hatano-Nelson model
We consider the thermodynamics of a non-Hermitian quantum gas with modes m and energies m . One can use a biorthogonal basis to compute the partition function of these systems, see Methods, where we introduce the subscript B for bosons and F for fermions, and ζ m = m − µ. From Z, one can obtain all the thermodynamic quantities, in analogy to what is done in a Hermitian system. Notice that these expressions have the same functional form as the ones for Hermitian gases [29]. However, there are important differences. The first is that the spectrum of non-Hermitian systems depends on the boundary conditions. This implies that the thermodynamic potentials will also change for different boundary conditions. In particular, the system with OBC is unstable for large system sizes, such that it can only be analyzed for small sizes. Nevertheless, the thermodynamic limit of this system can be achieved using the surrogate Hamiltonian [16], that has the same bulk spectrum. Furthermore, responses of the system, encoded in generalized forces (derivatives of the grand potential Ω with respect to the perturbation) are proportional to a biorthogonal correlation function [18]. More importantly, although the partition function is real in the presence of pseudo-Hermitian or parity-time (PT) symmetry [18] (see discussion in Methods), these complex energies change quantitatively the behavior of the thermodynamic potentials. In special, systems with Imζ m /Reζ m > 1 will show oscillations in β, together with a monotonic behavior.
To investigate this matter, we consider one paradigmatic non-Hermitian system, the Hatano-Nelson model [19,28] where M is the lattice size and c j (c † j ) annihilates (creates) a particle at site j. This is a simple hopping model with a reciprocal part (proportional to t) and a nonreciprocal part (proportional to Γ). The behavior of the spectrum depends on Γ/t and the boundary conditions, as discussed in detail in the Supplementary Material. In particular, for |Γ/t| > 1 and OBC, the system presents a completely imaginary spectrum, thus realizing an ideal platform to observe the oscillations in β. The results for several thermodynamic potentials are shown in Fig. 1. Figure 1: Thermodynamic potentials as a function of β for the Hatano-Nelson model with t = 0.1 for bosonic (red), fermionic (blue) and classical (green) systems. We set µ = −10 −3 , except for the bosonic system with PBC, where we use µ = −2t−10 −3 , such that the chemical potential is below the real part of the bands. For the periodic system we use Γ = 1, whereas for the open and surrogate systems we use Γ = √ 1 + t 2 , such that the typical energy scale ∆ associated with the Hatano-Nelson model is equal to 1 in all cases. The intensive grand potential f is shown in (a) for PBC, (b) surrogate, and (c) OBC. The intensive internal energy u is shown in (d) for PBC, (e) surrogate, and (f) OBC. The intensive entropy s is shown in (g) for PBC, (h) surrogate, and (i) OBC. The pronunced peaks that appear occur for 2∆ β = (2n + 1)π for fermions (shown with dashed blue lines in all figures) and 2∆ β = (2n)π for bosons (shown with dashed red lines in all figures). We use 1000 k-points for the periodic systems (equivalent to a lattice with 1000 sites) and a lattice with 20 sites for the open boundary system.
We consider the grand potential Ω, internal energy U and entropy S for the bosonic, fermionic, and classical realizations of the Hatano-Nelson model, see definitions in Methods. Some care should be taken with the values chosen for the chemical potential, as for bosons, µ < Re m for all modes [29]. Although the PT-broken phase for OBC displays purely imaginary energies, and hence this is not a problem, for bosons with PBC we need to choose µ < −2 |t|. Consequently, only energies at the minimum of the band will show oscillatory behavior, which will decay for increasing β. Nevertheless, for the other boundary conditions and for fermionic systems, the oscillations in β should be present, and µ will simply regularize the divergences of the thermodynamic potentials. The results are shown in Fig. 1 where n ∈ Z. There is, however, a noteworthy difference for the thermodynamic potentials for fermions with PBC and for bosons and fermions with OBC and the surrogate Hamiltonian. For PBC [Figs. 1 (a), (d) and (g)], only the large scale peaks are present. This is clearly seen when inspecting the internal energy, Fig. 1 (d), or the entropy, Fig. 1  , as they are both related to derivatives of f . The difference between both behaviors could have been anticipated because the spectrum for PBC has a non-zero real part for every Γ (see Supplementary Material), whereas the OBC system has a purely imaginary spectrum (see Supplementary Material) for |Γ| > |t|. Then, one can expect that more oscillations will be present.
These oscillations in β are very interesting because u is usually inversely proportional to β, such that an increasing temperature (decreasing β) would lead to a larger internal energy. Although there is a general decay of u for large β, which goes to zero in the limit of β → ∞ (T → 0), a small variation of β can lead to a large variation in the internal energy, a situation typical of the one occurring in the vicinity of a critical point [29].
The values of the peaks given by Eqs. (3) are the values that lead to the zeros (poles) of the fermionic (bosonic) partition function in Eq. (13), such that they describe a phase transition in the theory of Yang-Lee/Fisher zeros [27,[29][30][31][32][33][34][35]. Later, we will elaborate further on this special kind of phase transition.
The observation of oscillations in the thermodynamic quantities is not completely new. Intriguingly, features like those were also observed, although not discussed in these terms, in an analysis of a Wigner-Weyl representation of a non-reciprocal, therefore non-Hermitian, classical system [36]. Similarly, the Loschmidt overlap, which is the analogous of the free energy for evolution in real time, shows similar features for Hermitian dynamical phase transitions [37]. More interestingly, these oscillations are also a signature of the iTC phase proposed by Wilczek in his seminal paper on time crystals [20]. Those were studied in terms of dissipative systems in Ref. [38]. As we will show now, the conditions set by Eqs. (3) are precisely the ones that define the iTC phase.

Connection to imaginary time crystals
To investigate this matter, we express now the partition function as a path integral over coherent states in imaginary (or Euclidean) time τ [1,39,40]. This is done using the Trotter decomposition, where we introduce an (over) complete set of coherent states at every interval in imaginary time [39]. One can build such states using coherent states of a Hermitian operator, such as position, even for a non-Hermitian Hamiltonian. For completeness, we show in the Supplementary Material how to do second quantization and build coherent states for a pseudo-Hermitian operator. These states are parametrized by the fields Ψ and Ψ † and one can write the partition function as a path integral over them, where the Euclidean action S E is given (for a quantum gas) by with H the Hamiltonian matrix, such that the Hamiltonian density is Ψ † HΨ. These fields can depend on some continuum index, such as position or momentum, and can have also a tensorial structure accounting for spin or other inner degrees of freedom. Integration/summation/contraction over such degrees of freedom is implied. In addition, the components of Ψ will be Grassmann variables if it describes a fermionic field. Because for quantum gases the action is quadratic in the fields, one can exactly integrate this theory and the thermodynamic behavior is fully determined by ∂ τ − µ + H or, equivalently, the inverse of the Green's function G. As a matter of fact, [1] Besides determining the thermodynamic behavior of the system, G describes the response of the system to external perturbations. Wick theorem states [1] that any correlation function will be proportional to products of two-point functions, which are given by G. Therefore, G dictates the whole behavior of a quantum gas.
In particular, it is useful to consider the Fourier transform (FT) of G in τ . Due to the periodic conditions in τ , the frequencies will be discrete, being even (odd) multiples of π/β for bosons (fermions) [1,39]. Therefore, we have the Matsubara frequencies ω n ω n = π β n M = π β 2n, bosons, (2n + 1) , fermions, with n M the Matsubara mode. The FT of G, G n , is then given by where in the second equality we introduced twice the resolution of the identity m φ R m φ L † m =1 and used the biorthogonality of the eigenvectors φ R/L m of H with eigenstate m . Notice that H is an operator in position/momentum and a matrix on the inner degrees of freedom of the fields, but it is not an operator in the Hilbert space. As such, its eigenvectors φ R,L m (which are basically the wavefunctions) are just vectors of functions of position/momentum. We see Figure 2: Green's function and its FT for the fermionic Hatano-Nelson model at resonance with n M = ±11. Results are for PBC, surrogate Hamiltonian, and OBC for µ = −10 −3 , t = 0.1 and β = 1. We choose Γ such that the resonance condition is met (see main text) for all these boundary conditions. A different choice of β will only lead to a rescaling of G. We start by showing | G n | as a function of the Matsubara mode n M and k (distance to the left edge) for the periodic (open) system. A logscale is used for the periodic systems, such that the features besides the peaks of these functions are visible. This function is shown in (a) for the PBC system, in (b) for the surrogate Hamiltonian, and in (c) for the OBC system. Next, we show the real part of G(r, τ ) as a function of imaginary time τ and distance r, where the spatial dependence is either computed using a FT in momentum, in the case of the periodic systems, or is computed from the wavefunctions, for the open system. The results are shown in (d) for PBC, in (e) for the surrogate Hamiltonian, and in (f) for the OBC system. The results for the imaginary part are similar. We use 100 k-points for the periodic system and a lattice with 20 sites for the open system. then that the behavior of G n , and consequently G, is determined by ζ m . In particular, the poles of this function are given by which can only be satisfied (for real T or β) when Reζ m = 0 and Imζ m = ω n = n M π k B T . Notice also that ζ/(k B T ) should be an integer multiple of π. Thus, this is possible only in the low (or intermediate) temperature limit. When such conditions are met, both G B (τ ) and G F (τ ) take the form where ω n and ζ m almost satisfy Eq. 9. A small deviation from this condition, encoded in a very small but finite Reζ m , is necessary to prevent the function to diverge at the resonance. These systems then exhibit oscillatory behavior in imaginary time! Moreover, Eqs. It is remarkable that both G and G n show a direct signature of the iTC phase. We present the results for this function close to the resonance for the fermion mode n M = ±11 and µ = −10 −3 in Fig. 2. We plot | G F n | in Figs , with the difference that now they follow a sin(ka) function, changing the sign relation between n M and k. The OBC does not allow for an analysis in momentum space, and we must take into account effects that are not present for the periodic systems. One is the NHSE, which localizes φ R (φ L ) in the right (left) edge of the system. The other is the modification of the spectrum due to the small size of the system, making the bandwidth equal to 2  Fig. 2 (e)], similar patterns are visible but reversed, due to the fact that the poles in Fig. 2 (b) occur for opposite n M and k. These features are more visible for OBC [ Fig. 2 (f)], but with the complications commented before on the discussion for G n . Although there is a spatial decay due to the NHSE, the periodicity of this system does not change.
If Γ/(k B T ) is large, but not resonating, there will be much less intense poles at |n M | < 2Γ/π. This is clearly seen for PBC in Fig. 3 (a), where the resonances occur for |n M | = 1, 3, 5, but with a large spreading in momentum. Similar features are seen for the surrogate in Fig. 3 (b). Such peaks are shown for OBC in Fig. 3 (c), but with a stronger peak at n M = 3. Interestingly, the peaks in real space show different spatial periodicity, reflecting the fact that the peaks for different n M occur also for different k. This will influence the properties of G(r, τ ). The absence of a sharp peak leads to an incoherent behavior, Fig. 3  With these results, we understand which are the conditions for the occurrence of oscillations in imaginary time and real space. On resonance, Fig. 2, there is ordering in both imaginary time and real space. For low temperatures and off resonance, Fig. 3, many peaks will be present simultaneously, blurring the oscillations in τ . For high temperatures, Fig. 3, the oscillations are lost. The limit of T → 0 (β → ∞) in Eqs. (25) and (28) implies that these conditions will be satisfied for k = ±π/2a and all values of n M , such that the oscillations in τ will not be present. The behavior for the bosonic system is more intricate and is studied in the Supplementary Material. Now, we discuss the meaning of these phases in more depth. For Γ = 10k B T (βΓ = 10), the results for | G n | are shown in (a) for the PBC system, in (b) for the surrogate Hamiltonian, and in (c) for the OBC system. The real part of G(r, τ ) is shown in (d) for PBC, in (e) for the surrogate Hamiltonian, and in (f) for the OBC system. For Γ = 0.2k B T (βΓ = 0.2), the results for | G n | are shown in (g) for the PBC system, in (h) for the surrogate Hamiltonian, and in (i) for the OBC system. The real part of G(r, τ ) is shown in (j) for PBC, in (k) for the surrogate Hamiltonian, and in (l) for the OBC system. We use 100 k-points for the periodic system and a lattice with 20 sites for the open system.

Discussion
Time crystals are phases where the continuous time translation symmetry is reduced to a discrete translation symmetry, in analogy to what happens for spatial translation invariance in a crystal [20,41]. Even though the originally proposed model has some issues [42,43], it generated a substantial offspring [44][45][46]. From those, we highlight the recent proposal of time glasses [47], time quasicrystals [48] and dissipative time crystals [49][50][51][52], which are akin to the iTC phase.
The later was briefly conjectured at the end of Ref. [20], where the similarity between imaginary time and spatial dimensions in the Euclidean action are discussed. In analogy to the spatial variables, τ could also have preferred periods. In our case, Eq. (9) sets this period to be T = 2π/ω n = 2 β/(n M ) when n M is in resonance with a value of ζ. As the imaginary-time box has length β, for a resonance of Eq. (9), exactly n M /2 oscillations will be in the imaginary time interval, analogously to a standing wave in space. This is precisely what is seen in Figs The authors of Ref. [38] studied this phase looking at a bosonic system coupled to a bath. They had a non-local action in imaginary time, which corresponds to a non-Markovian evolution of the system and leads to an (imaginary) time-dependent Hamiltonian. Nevertheless, the results found there are similar to ours. Their model presents a charge-density wave order and the order parameter related to this phase shows oscillations in both β and τ . Eqs. (9) reveals that the peak for a specific Matsubara mode is also the peak for a specific k, see Methods. This will set a periodicity in space which, as seen from our results, survives also for a small OBC system, where momentum is no longer a good quantum number.
The existence of these oscillations in space is not incidental. Eq. (9) settles the condition for the presence of disorder lines for free systems [27] in the general theory of the Yang-Lee zeros. These phases were first obtained by Stephenson when studying the classical Ising model in a triangular lattice [21]. The correlation function of the order parameter has an oscillatory part, together with the exponential decay, typical of critical systems. The characteristic modulation length follows scaling laws [21][22][23][24][25][26] and is related to the presence of zeros of the partition function in the complex-parameter space. In the Euclidean action, τ is on the same footing as the spatial variables, so, intuitively, the presence of oscillations in imaginary time should not be surprising. The connection with disorder points also reveals that G n , for the value of n that presents a resonance, is a natural order parameter. The char-acterization of disorder points is done by using the Fourier transform of the correlation function [25,26] and is analogous to the definition of the order parameter in charge density waves [53].
However, imaginary time is distinct from real space in two aspects. First, bosons and fermions have different boundary conditions in τ , leading to different Matsubara frequencies for each of them. Second, the imaginary-time interval is set by temperature, being proportional to β, and τ is conjugated to energy. Therefore, the evolution in imaginary time, given by exp(−Hτ ), can be seen as a weighted projection on energy states [54]. For τ → ∞ (β → ∞), the system is projected to the ground state. The presence of a resonating condition implies that thare is a favored (imaginary) energy scale, which defines the period of oscillations in β. Due to pseudo-Hermiticity, the energies of the system come in complex conjugated pairs. Because the imaginary part of the energy is usually associated with dissipation, this can be interpreted as a kind of steady state, whereas the loss of probability in one mode is compensated by a gain in another. The fact that pseudo-Hermitian systems can have a biorthogonal unitary evolution [18] supports this view. Hence, such oscillations in temperature can be interpreted as a kind of resonating steady-state between system and reservoir, with the energy received/lost exhibiting peaks when the coupling of the system to the reservoir, given by Γ in this case, resonates with a Matsubara frequency, which is a typical frequency associated to the thermal state.
In this work, we studied the thermodynamics of non-Hermitian quantum gases for finite temperatures, with special focus on the Hatano-Nelson model. The model was chosen because it exhibits phases with purely imaginary single-particle energies for both PBC and OBC. The presence of such modes lead to the iTC phase conjectured by Wilczek in Ref. [20]. The existence of this phase is revealed by both, the oscillation in the thermodynamic potentials as a function of β, and by the oscillations of the Green's function as a function of imaginary time. There is also an order in real space because the iTC is a disorder point.
The iTC phase is interpreted as a resonance occurring precisely when the energy scale associated with the Matsubara frequencies matches the bandwidth of the imaginary part of the system, proportional to the coupling to the reservoir. Under this condition, the Matsubara frequency, usually interpreted as a mathematical tool, is manifested both in the thermodynamic and single-particle properties, and becomes measurable.
The results obtained here for the Hatano-Nelson model should be present also in more complicated non-Hermitian models, as long as they exhibit modes with purely imaginary energies. Different models will lead to different crystalline structures in imaginary time. The presence of multiple orbitals 13 and spin might unveil more intriguing phases. Studying these phases from the perspective of quantum heat engines and temperature dependent energy levels [55] may also reveal new and interesting heat phenomena.
An important remark is that the fact that the system is non-Hermitian ultimately comes from the interaction of the system with a reservoir. Therefore, the introduction of imaginary terms in the energies comes from a (zero frequency) imaginary part of a self-energy, similarly to what was considered in Refs. [38,56]. A frequency-dependent self-energy can lead to more intricate resonance conditions. A natural question concerning a more generic interacting system is the impact of many-body effects. In some discrete Floquet time-crystal, the effect of many-particle interactions is to stabilize the pre-thermal state, in which the time-crystalline behavior is observed, making it more robust [57]. Thermodynamic potentials are also multiparticle quantities and the effect of interaction can lead to combination of Matsubara frequencies in their observed behavior. Therefore, a further investigation of many-body effects in this description constitutes a promising topic for further research.
As the Hatano-Nelson model can be engineered in many platforms, the iTC phase can be observed experimentally. The detection of these effects require their realization in a quantum platform. The measurement of a thermodynamic quantity will give direct demonstration of the iTC phase. Measurements of correlation functions for different times (further analytically continued to imaginary time) should also yield an indication of this phase. Moreover, evolution in imaginary time is related to response to quantum quenches [54]. In this way, the response of the system to such a quench can also carry information on the periodicity in imaginary time.

Thermodynamic potentials of non-Hermitian quantum gases
For a non-Hermitian quantum gas, the energies may become complex and there are two kinds of eigenstates that label a microstate, |{n m } R and |{n m } L , which are eigenstates of H and H † , respectively, for each mode m with energy m ( * m ) and eigenstate |m R/L . The right eigenstates of H alone do not form a complete set. However, the left and right eigenstates together do form a complete basis [2,4,18]. Using these states as a biorthogonal basis, one can compute the grand partition function and obtain the usual result for a quantum gas whereN is the number operator, µ is the chemical potential (we will consider it to be real), and we define for convenience ζ m ≡ m − µ.
The partition function is, in general, complex, but it can be real in some special cases. In the presence of a pseudo-Hermitian symmetry, H and H † are related by a similarity transformation H † = gHg −1 , g = g † , such that their energies come in complex-conjugated pairs and their spectra are identical. As such, Z (computed from H) and Z * (computed from H † ) are the same, and consequently Z is real [18]. In addition, parity-time (PT) symmetry makes the energies to be real, and in this case Z is trivially real, as it is a sum of real numbers. In the following, it will become clearer that the most interesting effects for non-Hermitian gases will occur in the PT-symmetry broken phases, when the energies have a finite imaginary part, but the calculations presented here hold also if this symmetry is preserved.
For bosons, n m ∈ N, whereas n m ∈ {0, 1} for fermions. For bosons, we need to assume that Re ≥ µ (Re ζ m ≥ 0), |exp [−β n m ζ m ]| ≤ 1, otherwise this Laurent series diverges. The sum yields the familiar [29,58] partition functions for bosons and fermions, From Z, one can obtain all the thermodynamic quantities. Some thermodynamic potentials are particularly interesting to consider and are discussed in the Methods. The first is the average occupancy of each level: which leads to the Bose-Einstein and Fermi-Dirac distributions, Other kinds of potentials that are interesting are the ones related to the thermal behavior of the system. These are the grand canonical potential the internal energy U and the entropy S = −∂F/∂T = k B β 2 ∂F/∂β, which reads (18) where N = m N m is the total number of particles. Notice that this just follows from the thermodynamic definition of F [29]. These systems have a classical behavior in the limit |exp(−βµ)| 1 [29], when the length scale of thermal fluctuations (proportional to 1/T ) are smaller than the quantum ones (proportional to µ), which is related to high temperatures or small densities. In this situation, both functions reduce to the Boltzmann distribution The thermodynamic potentials then assume the form

Poles of the Hatano-Nelson model
For the system with OBC, one needs to numerically diagonalize H, to obtain the spectrum and the wavefunctions. As the Hatano-Nelson model has only one site per unit cell, for the periodic and surrogate Hamiltonians H will be just (k) and surr (k), respectively, which are numbers, and consequently do not have eigenvectors. The PBC system has the spectrum where a is the lattice parameter. If we just replace this (k) − µ in Eq. (9), we obtain the condition for resonance which is satisfied for Bosonic systems should satisfy µ ≤ −2 |t|. Therefore, for them the only possible solution for the above equations is k = 0 and n M = 0, which occurs when µ = −2t. Conversely, for fermions µ is not restrict and if we choose µ = 0, the resonance condition simplifies to k = ±π/(2a) and Γ = ∓n M (π/2)k B T . For OBC, a simple analysis is not really feasible. Hence, we turn to the surrogate Hamiltonian, as it has the same bulk spectrum. Using the band dispersion the condition for resonance becomes and we find the solutions In the PT-broken phase, |Γ| > |t|, and the first of Eqs. (28) can be rewritten as such that there will be poles at all values of n M that are smaller than the ratio 2 √ Γ 2 − t 2 / (πk B T ). For a finite system, however, k is of the form k = nπ/(M a), −M < n ≤ M , such that some of these modes can be present only for an infinite lattice, where k can take any value between −π/a and π/a.

Supplementary information
The Supplemental Methods contains details about the Hatano-Nelson model, discusses second quantization and path integrals for non-Hermitian operators and show the results for the Green's functions for bosonic systems.

Data availability
All data analyzed during the current study is available on reasonable request.  In this Section, we discuss in detail the Hatano-Nelson model. It is one of the earliest examples of non-Hermitian models [28] and is the simplest one to realize non-Hermitian topological phases [3]. In second quantization, its Hamiltonian reads

Author contributions
where M is the lattice size and c j (c † j ) annihilates (creates) a particle at site j. This is a simple hopping model in 1D with a reciprocal part (proportional to t) and a nonreciprocal part (proportional to Γ).
For non-Hermitian systems, the spectrum can be different for different boundary conditions. This is also the case for the Hatano-Nelson model. This occurs basically because of the breaking of Bloch theorem, as the wavefunctions are localized due to the non-Hermitian skin effect. One can then recover the bulk spectrum by doing an analytical continuation of the momentum, i. e., introducing an imaginary part that accounts for this localization. The Hamiltonian obtained then, the so called surrogate Hamiltonian [7,9,19], describes very well the bulk bands obtained for OBC.
In the case of the Hatano-Nelson model, the band dispersion for PBC is where a is the lattice parameter, while for the surrogate Hamiltonian it is (see Section 5.1.2) In Supplementary Figure 4, we show the spectrum of this system for both boundary conditions. For PBC, the real part of the spectrum is exactly the spectrum of a 1D tight-binding model, with a continuum of bands going from −2t to 2t. The system with OBC displays a more interesting behavior. For |Γ/t| < 1, it resembles also the spectrum of the simple hopping model, but with a bandwidth of 4 √ Γ 2 − t 2 instead of 4t. More surprisingly, for |Γ/t| > 1, the real part of the energy vanishes. Thus, for the system with OBC, the non-reciprocity caused by the coupling to the reservoir changes significantly even the real part of the spectrum.
The imaginary part of the spectrum is however the most interesting one in a non-Hermitian system. For the PBC system, there is again a continuum of bands, but now following sin(ka) instead of a cos(ka) and with a bandwidth directly proportional to Γ. Once again, the system with OBC shows more unconventional features, as there are finite imaginary values of the energy only for |Γ/t| > 1, when the energy is completely imaginary! The bandwidth is proportional to √ Γ 2 − t 2 , growing almost linearly with Γ for large Γ. Inspection of the spectrum shows an accumulation of modes at Γ = ±t, signaling the non-Bloch band collapse [9,16,[59][60][61]. At this point, the imaginary part of the momentum diverges (see Section 5.1.2), parity-time (PT) symmetry is broken, and the system is only pseudo-Hermitian, with the energy coming in complex conjugated pairs. For its simplicity and because it can show arbitrarily large imaginary parts in its spectrum (for both boundary conditions), the Hatano-Nelson is our model of choice to showcase the unique thermal features of non-Hermitian systems.

Surrogate Hamiltonian
From the second quantized version of the Hatano-Nelson model one can readily obtain this Hamiltonian in first quantization where |j denotes a one-particle state localized on site j. Then, the (time-independent) Schrödinger equation for a eigenstate |ψ R with energy takes the form where we used that the wavefunction on site j is j|ψ R = φ R (j) and j|l = δ j,l .
One can now assume that the wavefunction has a non-Bloch form [7][8][9] φ R (l) = e i[k( )+iκ( )](l−j)a φ R (j), with an imaginary part κ of the momentum, which will take into account the localization of the wavefunction. Then, one can convert the Schrödinger equation in an equation for k, κ and , Assuming that κ does not depend on , but only on the parameters of the model, this equation simplifies to the dispersion relation of the surrogate Hamiltonian 25 We can choose a specific value of to obtain κ. Choosing = 0 Substituting then the expression of exp(−κa) in Eq. (37), one obtains If t > 0 and Γ > t, and surr = −2 |Γ 2 − t 2 | cos(ka). For negative t, a similar argument holds, but with a difference in sign in the last equality. We can then write the expression for surr (k) as which is precisely Eq.

Biorthogonal second quantization and coherent states
Let us define the eigenstate |m R as being given by a creation operator R † m acting on the vacuum state |0 R , and a similar relation holds between |m L and a creation operator L † m , These eigenstates must satisfy the biorthogonal relations. A very natural way to do it, is to impose (anti)comutation relations [62,63] between R † and L where we use that L 0|0 R = 1, that L m annihilates |0 R , L m |0 R = 0, and that the commutator L m , R † such that and Similar relations hold between R and L † , In addition, we have that such that the action of R † and L (L † and R) on many-body right (left) states is the same as regular creation and annihilation operators [39,64]. In particular, the eigenvalues of the right (left) number operator N R m (N L m ) of a mode m, defined as R † m L m (L † m R m ), are the natural numbers. For fermions, only 0 and 1.
Using these properties, we can define coherent states [39] of a given mode with the upper sign for bosons and the lower one for fermions, which are eigenstates of L m and R m , Note that for fermions, ψ m is a Grassmann variable. We can then define a coherent state for all modes where Ψ = ψ 0 ψ 1 . . . is a vector of the ψ m parameters. The coherent states are not biorthogonal but they do form a supercomplete set and we can use them as a basis to write any operator. In particular, where Ψ is composed of Grassmann variables for fermions and is a vector of complex number for bosons [39].

Exact expresssions for G(r, τ ) for a free system
The Fourier transform of the Green's function determines G(r, τ ). One can perform the Matsubara sums exactly (done in Mathematica [65]) in G n to obtain G B and G F to denote bosonic (n M even) and fermionic (n M odd) functions. Notice that in the presence of poles these functions also reduce to the form as in the main text.

iTC for bosons
In this Section, we show the Green's function for the bosons. The fact that such systems have peaks for n M = 0 make that purely spatial oscillations appear together with the plane waves in imaginary time, seen for the fermionic systems. We analyze G B n and G(τ, r) B for Γ: (i) resonant, Supplementary The G n for PBC at a resonance (Supplementary Figure 6), presents a sequence of small peaks together with a very bright peak at n M = 0 and k = 0. For the surrogate, there will be peaks at n M = ±10 and k = ±π/2, together with the ones at n M = 0 for k = 0 or k = π. For OBC, there are only peaks at n M = ±10, without any signature at n M = 0. This is a finite-size effect, as small systems have a gap of order O (1/M ). As the case for the fermions, the pole structure of G n explains the features in G(r, τ ), with the PBC system having an amorphous behavior, the surrogate showing an interference pattern, and the OBC system showing behavior typical of a wavefront.
For large, non-resonating Γ (Supplementary Figure 7), the system presents pronounced peaks only at n M = 0 for PBC and surrogate, and at n m = ±4 for the OBC. The fact that the peak at n M = 0 occurs for only k = 0 for PBC makes that G(r, τ ) has features of an amorphous phase. Conversely, for the surrogate, there is a peak also at n M = 0 and k = π/a, so there will be oscillations on space with λ = 2π/k = 2a. The presence of the peaks at n M = ±4, favor oscillations with T = 0.5 β, although they are not ordered because there is a significant weight on other modes.
For small, non-resonating Γ (Supplementary Figure 8), there are pronounced peaks only for n M = 0. For the PBC and the surrogate, the weight of the poles is distributed among different values of k. For OBC, there are oscillations in space but a very quick exponential decay of the correlation function. There is modulation in space for PBC and surrogate, although no predominant wavelength is observed. For OBC, there is only decay in space.
Supplementary Figure 6: Green's function and its FT for the bosonic Hatano-Nelson model at resonance with n M = ±10. Results are for PBC, for the surrogate Hamiltonian, and for OBC, for t = 0.1, µ = −2t − 10 −3 (PBC) or µ = −10 −3 (OBC and surrogate), and β = 1. We choose Γ such that the resonance condition is met (see main text) for all these boundary conditions. A different choice of β will only lead to a rescaling of G.We start by showing | G n | as a function of the Matsubara mode n M and k (distance to the left edge) for the periodic (open) system. A logscale is used for the periodic systems, such that the features besides the peaks of these functions are visible. This function is shown in (a) for the PBC system, in (b) for the surrogate Hamiltonian, and in (c) for the OBC system. Next, we show the real part of G(r, τ ) as a function of imaginary time τ and distance r, where the spatial dependence is either computed using a FT in momentum, in the case of the periodic systems, or is computed from the wavefunctions, for the open system. The results are shown in (d) for PBC, in (e) for the surrogate Hamiltonian, and in (f) for the OBC system. The results for the imaginary part are similar. We use 100 k-points for the periodic system and a lattice with 20 sites for the open system. Supplementary Figure 7: Green's function and its FT for the bosonic Hatano-Nelson model for a large, but non-resonant Γ. Results are for PBC, for the surrogate Hamiltonian, and for OBC, for t = 0.1, µ = −2t − 10 −3 (PBC) or µ = −10 −3 (OBC and surrogate), Γ = 10k B T (βΓ = 10) and β = 1. A different choice of β will only lead to a rescaling of G.We start by showing | G n | as a function of the Matsubara mode n M and k (distance to the left edge) for the periodic (open) system. A logscale is used for the periodic systems, such that the features besides the peaks of these functions are visible. This function is shown in (a) for the PBC system, in (b) for the surrogate Hamiltonian, and in (c) for the OBC system. Next, we show the real part of G(r, τ ) as a function of imaginary time τ and distance r, where the spatial dependence is either computed using a FT in momentum, in the case of the periodic systems, or is computed from the wavefunctions, for the open system. The results are shown in (d) for PBC, in (e) for the surrogate Hamiltonian, and in (f) for the OBC system. The results for the imaginary part are similar. We use 100 k-points for the periodic system and a lattice with 20 sites for the open system. Supplementary Figure 8: Green's function and its FT for the bosonic Hatano-Nelson model for small Γ. Results are for PBC, for the surrogate Hamiltonian, and for OBC, for t = 0.1, µ = −2t − 10 −3 (PBC) or µ = −10 −3 (OBC and surrogate), Γ = 0.2k B T (βΓ = 0.2) and β = 1. A different choice of β will only lead to a rescaling of G.We start by showing | G n | as a function of the Matsubara mode n M and k (distance to the left edge) for the periodic (open) system. A logscale is used for the periodic systems, such that the features besides the peaks of these functions are visible. This function is shown in (a) for the PBC system, in (b) for the surrogate Hamiltonian, and in (c) for the OBC system. Next, we show the real part of G(r, τ ) as a function of imaginary time τ and distance r, where the spatial dependence is either computed using a FT in momentum, in the case of the periodic systems, or is computed from the wavefunctions, for the open system. The results are shown in (d) for PBC, in (e) for the surrogate Hamiltonian, and in (f) for the OBC system. The results for the imaginary part are similar. We use 100 k-points for the periodic system and a lattice with 20 sites for the open system.