A new type of non-Hermitian phase transition in open systems far from thermal equilibrium

We demonstrate a new type of non-Hermitian phase transition in open systems far from thermal equilibrium, which can have place in the absence of an exceptional point. This transition takes place in coupled systems interacting with reservoirs at different temperatures. We show that the spectrum of energy flow through the system caused by the temperature gradient is determined by the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi^{4}$$\end{document}φ4-potential. Meanwhile, the frequency of the maximum in the spectrum plays the role of the order parameter. The phase transition manifests itself in the frequency splitting of the spectrum of energy flow at a critical point, the value of which is determined by the relaxation rates and the coupling strength. Near the critical point, fluctuations of the order parameter diverge according to a power law with the critical exponent that depends only on the ratio of reservoirs temperatures. The phase transition at the critical point has the non-equilibrium nature and leads to the change in the energy flow between the reservoirs. Our results pave the way to manipulate the heat energy transfer in the coupled out-of-equilibrium systems.

www.nature.com/scientificreports/ the spectrum of energy flow between the reservoirs. We show that this spectrum is determined by an analog of the ϕ 4 -potential, and that it exhibits frequency splitting at a critical point (CP) above which there are two minima in the potential. Using an analogy with the theory of second-order phase transitions, we determine the frequency of the maximum in the spectrum as an order parameter for this non-Hermitian transition. We demonstrate that near the CP, fluctuations of the order parameter diverge according to a power law. The corresponding critical exponent depends only on the ratio of the temperatures of the reservoirs, and remains unchanged when the temperatures of all the reservoirs change by the same factor. This indicates the non-equilibrium nature of the phase transition at the CP. It is remarkable that although the energy flow between the reservoirs is proportional to the difference of their temperatures, the CP does not depend on the reservoir temperatures. Like an EP, the CP depends only on the relaxation rates of the system and the coupling strength between its subsystems. Our results open the way for observing non-Hermitian phase transitions in systems far from thermal equilibrium. The conditions for a CP in this transition differ from the conditions for an EP, thus making it possible to observe non-Hermitian phase transitions in systems without an EP, e.g., in a system of two coupled oscillators with the same relaxation rate but with different reservoirs temperatures. The predicted non-Hermitian phase transition leads to the change in the energy flow between the reservoirs that paves the way to control the heat energy transfer in nanoscale and can be used to develop phonon devices [36][37][38]40,41 .

Model
We consider two coupled oscillators interacting with their own reservoirs. The reservoirs are in thermodynamic equilibrium with temperatures T 1 and T 2 , respectively. This system can be implemented, for example, based on optically coupled nanomechanical resonators 42 . By eliminating the degrees of freedom of the reservoir using the Born-Markovian approximation 43,44 , we can obtain the equations for the oscillator slow amplitudes 43,44 (see also "Methods" section): Here, � a = (a 1 , a 2 ) T is a vector of the amplitudes of the first and second oscillators, respectively; is the coupling strength between the oscillators; γ 1,2 are the relaxation rates of the oscillators; and � ξ = ξ 1 ξ 2 T is a noise term that always appears together with relaxation terms 43,44 . In the Born-Markovian approximation, the noises are delta-correlated and are connected with the relaxation rates according to 43,45 where D = γ 1 T 1 0 0 γ 2 T 2 is a diffusion matrix (we put k B = 1 ). For brevity, we denote the matrix on the right- We calculate the spectrum of stationary fluctuations (i.e., at times t ≫ γ −1 1,2 ), of the system using the Wiener-Khinchin theorem 46,47 : where � a * (τ )� a T st is a matrix of two time correlators that is determined using regression theorem as follows (see "Methods" section for details): The averages � a * � a T st = � a * (τ = 0)� a T st obey the following condition 43 (see also the "Methods" section): It should be noted that the noise results in nonzero stationary values for the oscillator energies a * 1 a 1 st and a * 2 a 2 st as well as the value a * 1 a 2 st . The real part, Re a * 1 a 2 st , represents the interaction energy between the oscillators, while the imaginary part, Im a * 1 a 2 st , represents the energy flow from the first oscillator to the second. Using the Eqs. (3) and (4) we obtain the following expression for the spectrum 43 : Using the definitions of matrices M and D , we obtain: cr )ω 2 + ω 4 and � 2 cr = γ 2 1 + γ 2 2 /2 . The diagonal elements of Ŝ (ω) determine the spectra of the fluctuations of the first and second oscillators. The non-diagonal elements determine www.nature.com/scientificreports/ the spectrum of interaction between the oscillators. The real part determines the spectrum of fluctuations of interaction energy, while the imaginary part determines the spectrum of energy flow from one oscillator to another.

Spectra of oscillators
In the absence of noise, the stationary state of the system of Eq. (1) is zero. The dynamics of the system when it tends to the stationary state is determined only by the eigenvalues and eigenstates of matrix M , which have the form: and In this system, there is an EP at which the eigenvalues are equal to each other and the eigenstates coincide. The EP arises when � = � EP = |γ 2 − γ 1 |/2 , and separates the weak ( � < � EP ) and strong ( � > � EP ) coupling regimes 3,4 .
It is assumed that the transition to the strong coupling regime can be detected based on the splitting in the system spectrum 3,4 . However, splitting in the fluctuation spectra of the first and second oscillators occurs at coupling strengths of split 1 and split 2 , which differ from the coupling strength at the EP ( EP ) (Fig. 1) and from each other (the analytical formulas for split 1,2 are given in the "Methods" section). That is, the splitting in the spectra of the first and second oscillators occurs at different values of the coupling strength 42 . split 1,2 depend on the ratio of the reservoir temperatures ( T 2 /T 1 ) (Fig. 2). It should be noted that the splitting in the oscillator spectra can take place even when � < � EP , i.e., in the weak coupling regime (see the dashed red line in Fig. 2).  www.nature.com/scientificreports/ Thus, the splitting point in the system spectrum depends on the reservoir temperatures, and generally does not coincide with the EP. We note that at this stage, it is difficult to provide a direct analogy with the standard theory of second-order phase transition, due to absence of an order parameter and the divergence of its fluctuations near the transition point. As we will see below, this analogy can be established if we consider the spectrum of the energy flow.

Non-Hermitian phase transition in the spectrum of energy flow between the reservoirs
Energy flow between the reservoirs. As mentioned above, there is an energy flow through the system, i.e., the system is not in thermodynamic equilibrium. The spectrum of this energy flow is determined by the product of the imaginary parts of the non-diagonal elements of the matrix Ŝ (ω) and the coupling strength : cr )ω 2 + ω 4 determines the frequency spectrum of J(ω) , which splits when � 2 = � 2 cr = γ 2 1 + γ 2 2 /2 . Indeed, the quantity J(ω) represents the average energy flow between the reservoirs, and is always directed from a hotter to a colder reservoir (Eq. (9)). The amplitude of J(ω) depends on the difference of the reservoir temperatures. However, the form of the spectrum J(ω) does not depend on the reservoir temperatures, and is determined only by the denominator �(ω) . When � 2 < � 2 cr , there is a single maximum in the spectrum J(ω) at ω max = 0 . In the opposite case, when � 2 > � 2 cr , there are two maxima in the spectrum J(ω) at ω max = ± � 2 − � 2 cr . At 2 = 2 cr , splitting arises in the spectrum of the energy flow. Note that the integral of imaginary part of the non-diagonal element of the matrix Ŝ (ω) over all frequency range, �ImS 12 � = +∞ −∞ �ImS 12 (ω)�dω , depends non-monotonically on the coupling strength and the ratio of the relaxation rates, γ 2 /γ 1 . If γ 1 = γ 2 , ImS 12 is maximum when 2 = 2 cr , i.e., at the splitting point in the spectrum of energy flow (see Fig. 3a). This quantity reaches its absolute maximum under the condition γ 1 = γ 2 and 2 = 2 cr (see Fig. 3b), i.e., when � = γ 1 = γ 2 = � 0 cr . Note that cr = EP and, moreover, at the = cr , the eigenvalues and eigenstates of the system matrix M do not qualitatively changes. However, in the same way as EP , cr does not depend on the reservoir temperatures. The expression for �(ω) is an analog of the ϕ 4 -potential for order parameter 48 , while ω max is an analog of an order parameter at the second-order phase transition. As we will show below, the fluctuations of ω max diverge near cr .

Critical behavior near the non-Hermitian transition.
To establish the similarity between the non-Hermitian phase transition in a non-equilibrium system and the second-order phase transition, we study the fluctuation behavior of energy flow J(ω) near the critical coupling strength.
On average, the energy flow is directed from a hotter to a colder reservoir. The spectrum of the average energy flow, J(ω) , is determined by Eq. (9). To study the fluctuations of the energy flow, we simulate the system dynamics using Eq. (1) with noise. We calculate the stochastic time evolution of the energy flow over a finite time range t ∈ 0, t f , and then find the Fourier spectrum of the energy flow, J(ω) . For this purpose, we simulate the Eq. (1) using the Euler difference scheme that is stable near the stationary state 49,50 . The noise terms acting on the first and second oscillators are modeling as independent of each other random processes with independent increments over time. These noises are delta-correlated in time (white noise) and their average values are zero, see Eq. (2) and "Methods" section. The noise correlators are determined by diagonal elements of the diffusion www.nature.com/scientificreports/ matrix D , which are proportional to the product of the reservoir temperature and the relaxation rate of the corresponding oscillators (see Eq. (2), "Methods" section) 50 . By averaging over a large number of simulations of Eq. (1), we obtain spectra for J(ω) using Eq. (9) (see Fig. 4). In a given realization, the energy flow can fluctuate and the direction of flow may change. As a result, the spectrum of energy flow calculated from a single realization differs significantly from J(ω) (Fig. 5). However, the frequency distribution of the maxima in the spectrum J(ω) resembles the form of 1/�(ω) (cf. the solid blue and dashed red lines in Fig. 5). That is, �(ω) plays the role of a potential for the distribution of maxima in the spectrum J(ω) . Accordingly, the frequency of the maximum in the spectrum J(ω) can be considered as the order parameter of the non-Hermitian phase transition in a system far from equilibrium. In turn, the splitting point in the spectrum J(ω) can be referred as a critical point (CP) of the phase transition.
To describe the fluctuations of the order parameter, we find the frequency of the maximum in the spectrum J(ω) , i.e., ω max , for each of the simulations. We then average ω max and ω 2 max over a large number of simulations, and calculate the dispersion D(ω max ) = ω 2 max − �ω max � 2 for different values of the coupling strength and the reservoir temperatures T 1,2 (see Fig. 6a). We can see that the behavior of D(ω max ) changes qualitatively at the point of phase transition ( 2 = 2 cr ). Near the critical point (CP), the dispersion behaves as D(ω max ) ∼ |� − � cr | α . In accordance with the terminology used in the theory of second-order phase transitions 48 , we refer to α as the critical exponent. It can be shown that the dispersion D(ω max ) and the critical exponent α depend on the ratio of the reservoir temperatures ( T 2 /T 1 ) (see Fig. 6b), but do not depend on the absolute values of these temperatures (see "Methods" section).
Note that within a small neighborhood of the transition point, the power law is violated due to the finiteness of the system, a result that also agrees with the theory of second-order phase transitions 48 .

Conclusions
We demonstrate that in non-Hermitian systems far from equilibrium, a new type of non-Hermitian phase transition takes place. We consider a system of two coupled harmonic oscillators interacting with reservoirs at different temperatures. The difference of the reservoir temperatures moves the system away from thermal equilibrium, resulting in an energy flow through the system. We show that the spectrum of this energy flow is determined by an analog of the ϕ 4 -potential, and exhibits frequency splitting at a CP above which there are two minima in the potential. This spectral splitting at the CP can be associated with a non-equilibrium phase transition. The frequency of the maximum in the spectrum plays the role of an order parameter for this phase transition. We demonstrate that near the CP, fluctuations of the order parameter diverge according to a power law, and we  www.nature.com/scientificreports/ calculate the critical exponent. We show that this depends only on the ratio of the reservoir temperatures, which indicates the non-equilibrium nature of the phase transition. The non-Hermitian phase transition described here can take place in non-Hermitian systems without an EP. Moreover, this phase transition is not accompanied by qualitative changes in the eigenvalues and eigenstates of the system matrix. This opens the way for the study of non-Hermitian phase transitions in a new class of systems.

Methods
Derivation of the system equations. We consider the system of two coupled oscillators interacting with their reservoirs, which are described by infinite sets of oscillators. The Hamiltonian of this system is is a Hamiltonian of the system of two coupled oscillators in the rotating-wave approximation 46 ( = 1 ), where â 1,2 and â † 1,2 are the annihilation and creation operators of the first and second oscillators, obeying the boson commutation relations 46 . ω 0 is a frequency of individual oscillators. is a coupling strength between the oscillators. Ĥ R1 = , b † 2k 2 are the annihilation and creation operators of the oscillators in the first and second reservoirs. ω 1k 1 and ω 2k 2 are frequencies of these oscillators.
are Hamiltonians of the interaction between the first/second oscillators and the first/second reservoir, respectively. g 1k 1 and g 2k 2 are the coupling strength between the first/second oscillator and the k 1,2 -th oscillator in the first/second reservoirs.
Using the Heisenberg equation for operators 43,44 and then moving from the operators to averages, we obtain the following equations Here a 1,2 = â 1,2 , b 1k 1 = b 1k 1 , b 1k 2 = b 1k 2 are averages of the respective operators. Then, we formally integrate the last two equations: www.nature.com/scientificreports/ Substituting Eqs. (15) and (16) into Eqs. (11) and (12), respectively, we obtain Making a change of variables ã 1,2 = a 1,2 e iω 0 t , we obtain the equations for slow variables: Following the Born-Markov approximation 43,44 and using the Sokhotskii-Plemelj formula 51 , we calculate the integrals in the right parts of Eqs. (19) and (20). As a result, we obtain the equations with the relaxation and noise terms: Considering that the first and second reservoirs are in thermal equilibrium with temperatures T 1 and T 2 , respectively, we find the correlation properties of noise: where � ξ(t) = (ξ 1 (t), ξ 2 (t)) T is a vector of noise terms and D = γ 1 T 1 0 0 γ 2 T 2 is a diffusion matrix 43,44 . Thus, it is seen that in the derived equations, the relaxation and noise terms are related to each other by the fluctuationdissipative theorem 43,44 . Derivation of the expression for the system spectrum. Using the Wiener-Khinchin theorem 46,47 , we calculate the spectrum of a system We first calculate the correlator � a * (τ )� a T st in Eq. (24). To do this, we formally integrate the matrix Eq. (1): Then, using the quantum regression theorem 43 , we obtain: www.nature.com/scientificreports/ All eigenvalues of the matrix M have negative real parts. Thus, at times t ≫ γ −1 1,2 , the first term on right side in Eq. (26) tends to zero. Taking into account the expressions for the noise correlators (see Eq. (2)) and considering the case t ≫ γ −1 1,2 , we derive: When τ → +0 , Eq. (27) takes the form: Differentiating both sides of Eq. (28), we obtain: In the steady state, � a * (t)� a T (t) is determined by the following equation: The expression in Eq. (31) holds true when τ > 0 . To use the Wiener-Khinchin theorem 46,47 we need to calculate the correlator at τ < 0 . This correlator is calculated in the same way, as follows: Thus, the spectrum can be expressed as: By multiplying the spectrum matrix on the left and right by the respective matrices and taking into account Eqs. (31) and (32), we obtain: Finally, we obtain an expression for the spectrum matrix: where �(ω) = (γ 1 γ 2 + � 2 ) 2 − 2(� 2 − � 2 cr )ω 2 + ω 4 and � 2 cr = γ 2 1 + γ 2 2 /2. The diagonal elements of Ŝ (ω) define the spectra of the first and second oscillators, and the non-diagonal elements are complex quantities. The real parts of the non-diagonal elements are equal to each other, and define the spectrum of interaction between the oscillators, while the imaginary parts of the non-diagonal elements differ from each other in sign and define the spectra of the energy flow from the first oscillator to the second and vice versa.
Note that ω is the difference between the oscillation frequency and the frequency of a single oscillator ω 0 (we consider that the frequencies of oscillators are equal to each other).
Dependence of the dispersion of the order parameter on the ratio of the reservoir temperatures. Numerical simulation of Eq. (1) shows that the spectra of the oscillators, the spectra of the interaction and the energy flow between the oscillators, and the dispersion of the frequency of the maximum in the spectrum ( ω max ) of energy flow J(ω) depend on the ratio of the reservoir temperatures rather than their absolute values. To ascertain the mechanism of this dependence, we analyze the behavior of the system with noise, as follows: where � a = a 1 a 2 T is a vector of the amplitudes of the first and second oscillators, respectively; is the coupling strength between the oscillators; γ 1,2 are the relaxation rates of the oscillators; and � ξ = ξ 1 ξ 2 T is a noise term that obeys the following conditions In the stationary state, the amplitudes a 1,2 are nonzero only when at least one of the reservoir temperatures is nonzero. Without loss of generality, we assume that T 2 = 0 . In this case, we can make changes to the variables ã 1,2 = a 1,2 / √ T 2 and ξ 1,2 = ξ 1,2 / √ T 2 . In these new variables, the Eq. (37) can be rewritten as: The redefined noise terms obey the conditions where D = γ 1 T 1 /T 2 0 0 γ 2 is a redefined diffusion matrix. It is important to note that D , and consequently the amplitudes � a = ã 1ã2 T , depend only on the ratio of the reservoir temperatures. As a result, the spectrum matrix S (ω) calculated by Eq. (39) depends on the ratio of the reservoir temperatures (where the method of calculation is the same as for Eq. (37)). The spectrum matrix Ŝ (ω) is determined by S (ω) as (see the definition of Ŝ (ω) in Eq. (3)): Thus, the form of spectrum Ŝ (ω) depends only on the ratio of the reservoir temperatures, while the absolute value of the temperature T 2 is determined only by the amplitude of the spectrum. The same conclusion holds true for a spectrum calculated based on a single simulation of Eq. (37). As a result, the frequency of the maximum in the spectrum ( ω max ) calculated for this simulation also depends only on the ratio of the reservoir temperatures. Since the dispersion of the frequency D(ω max ) = ω 2 max − �ω max � 2 is calculated by averaging ω max and ω 2 max over a large number of simulations, then D(ω max ) depends only on the ratio of the reservoir temperatures. (36) S(ω) = 1/π �(ω) γ 1 T 1 (γ 2 2 + ω 2 ) + γ 2 T 2 � 2 i �(γ 2 T 2 (γ 1 + iω) − γ 1 T 1 (γ 2 − iω)) i� (γ 1 T 1 (γ 2 + iω) − γ 2 T 2 (γ 1 − iω)) γ 2 T 2 (γ 2 1 + ω 2 ) + γ 1 T 1 � 2 , (37) d dt a =M� a +� ξ(t).