Protection of Correlation-Induced Phase Instabilities by Exceptional Susceptibilities

At thermal equilibrium, we find that generalized susceptibilities encoding the static physical response properties of Hermitian many-electron systems possess inherent non-Hermitian (NH) matrix symmetries. This leads to the generic occurrence of exceptional points (EPs), i.e., NH spectral degeneracies, in the generalized susceptibilities of prototypical Fermi-Hubbard models, as a function of a single parameter such as chemical potential. We demonstrate that these EPs are necessary to promote correlation-induced thermodynamic instabilities, such as phase-separation occurring in the proximity of a Mott transition, to a topologically stable phenomenon.

Introduction.Topology has been introduced in physics to understand robustness.A pioneering achievement of this approach is the understanding of the quantum Hall effect [1], where the striking quantization of a transverse conductance has been explained in terms of topological properties of Bloch bands [2][3][4].Since then, topology has conquered a wide range of physical settings far beyond the band theory of solids [5][6][7][8][9][10][11][12][13][14][15].As a novel direction within this paradigm, here we explain the robustness of phase instabilities in Fermi-Hubbard models by revealing and studying unique topological properties of their generalized susceptibilities [16,17].These quantities, which encode information on the electronic fluctuations of a many-body system, can be linked to basic observables such as the uniform charge response (isothermal compressibility) χ q=0 [16,18].
Importantly, even for closed many-body systems described by a Hermitian Hamiltonian, such generalized susceptibilities naturally acquire a complex spectrum as matrices in Matsubara frequency space, and are thus subject to a non-Hermitian (NH) topological classification approach [15,19].There, we identify the emergence of inherent NH symmetries [20,21], requiring merely the physical assumptions of thermal equilibrium and Fermi statistics of the constituents.Under these ubiquitous circumstances, pairs of exceptional points (EPs), i.e., NH spectral degeneracies at which the generalized charge susceptibility matrix χ νν ′ c becomes non-diagonalizable [22][23][24][25][26][27][28][29], generically occur as a function of a single tuning parameter such as chemical potential (or filling fraction).These EPs are topologically protected by the aforementioned inherent NH symmetry and their splitting in parameter space (see Fig. 1

(b) for an illustration).
As a consequence, purely real eigenvalues λ I of χ νν ′ c occur in an extended parameter range in between the EPs and generically trigger divergences in the uniform charge response χ q=0 [30] that are robust against small parameter changes (cf.Fig. 1(a) for an illustration).These FIG. 1. Schematic illustration of the enhancement of the uniform charge response χq=0 (a) triggered by the eigenvalue λI reaching the condition λI = −1/t 2 (b).The chemical potential µ at which the divergence happens is located within a finite range, bounded by the EPs (red dots on (b)).If a divergence is found, no small perturbation can wash it out, as the intersection between the "rectified" λI curve (thin gray line) and the −1/t 2 limit can only be continuously shifted.
divergences in χ q=0 signal the propensity of the correlated system to undergo a thermodynamic phase separation between a compressible metallic and an almost incompressible "bad metal" phase, often occurring in the proximity of Mott metal-to-insulator transitions [30][31][32][33][34].This phase separation can be viewed, in many respects, as the electronic counterpart of the liquid gas transition for water molecules.
Non-Hermitian topology and EPs have been widely discussed in systems where the Bloch band structure has been augmented by dissipative terms of various physical origin, ranging from scattering rates of quasi-particles [35][36][37][38][39][40][41][42][43][44][45] to gain and loss in optical systems [46][47][48][49][50][51][52].There, the experimental visibility of NH signatures is oftentimes limited by the overall blurring arXiv:2307.00849v3[cond-mat.str-el]7 May 2024 introduced by imaginary damping terms in the dissipative time-evolution.By contrast, since the generalized susceptibilities at the heart of our present analysis do not relate to effective complex energy spectra, their NH topology has a direct impact on natural observables, independent of idealistic assumptions on temperature and without the need for complex multi-orbital models, respectively.
NH symmetries of generalized susceptibilities.The central objects of our analysis are four-point functions describing the propagation of a particle-hole pair (see Fig. 2(b)) in the setting of a time-independent Hamiltonian H at thermal equilibrium.These are expressed as matrices of two fermionic Matsubara frequencies ν and ν ′ , where ν (′) = (2n (′) + 1)π/β, n (′) ∈ Z, and β = 1/T is the inverse temperature.In the literature, such quantities are referred to as generalized susceptibilities, since they yield static physical response functions when summed over both fermionic frequencies [53].Specifically, we define them as (illustrated in Fig. 2 αi e −Hτ the Fourier transform of the (creation) annihilation operators [54], and G (2)νν ′ α1...α4 the two-particle Green's function.
α i refers to, in principle, any of the degrees of freedom of the model (momenta, spin, orbital, etc.).Some properties of χ νν ′ ph,α1...α4 have been already analyzed in Refs.[16,17,53,55,56].In this work, we are investigating the topological properties of the corresponding eigenvalue spectrum.
, where G(iν) refers to the oneparticle Green's function (here, we dropped the αi indices for clarity).Panel (c), schematic illustration of the centro-Hermitian matrix symmetry.
In the following, we consider the local generalized charge susceptibility ph,σσσ ′ σ ′ of a oneorbital model that satisfies the following relation: and is therefore a centroHermitian matrix.In addition, if the Hamiltonian possesses specific symmetries, these can impose even stricter matrix properties.For instance, for particle-hole symmetry (PHS) χ νν ′ c becomes real and has only real eigenvalues [17,56].
Minimal model for exceptional susceptibilities.To illustrate our general findings, we consider a 2 × 2 matrix χ 2×2 ph obeying the centroHermitian condition: where a, b, c, d ∈ R and is a complex vector.The a parameter can be safely disregarded, as it only amounts to a rigid eigenvalue shift.EPs are globally stable for a two-or higher-dimensional parameter space, because for the matrix to become nondiagonalizable, two conditions (v 2 R − v 2 I = 0, ⃗ v R • ⃗ v I = 0) have to be simultaneously satisfied.It is immediate to see that the centroHermitian property implies that the second is always fulfilled.It is then sufficient, for the exceptional points to manifest, that c 2 + d 2 − b 2 = 0, which implies that even in a one-dimensional space, EPs -if any are present -will be globally robust against any perturbation representable by a matrix of the form given in Eq. ( 4) [61].Crucially, no other perturbation can arise, because the centroHermitian condition does not originate from any further symmetry, but it is an intrinsic consequence of the time-independence of H and a defined quantum (here: Fermi-Dirac) statistics [62].The occurrence of these stable EPs generically mark the borders of regions with complex conjugate pairs on the one hand and real and distinct eigenvalues on the other hand.On the contrary, if we additionally impose PHS, χ 2×2 ph becomes purely real and symmetric (b, d = 0), which implies that the only solution of the two conditions for the EPs will be for ⃗ v R = ⃗ v I = 0.This is known in literature as a diabolic point, which is effectively concurrent with a Hermitian degeneracy [15] and, indeed, generally requires fine-tuning.Analytical study of the atomic limit.As the simplest physical platform to exemplify the spectral properties of χ νν ′ c , we now study the exactly solvable atomic limit of the Hubbard model (AL) where µ is the chemical potential, n σ = c † σ c σ the occupation of an electron with spin σ, and U the on-site Coulomb repulsion given in arbitrary units of energy.This model fulfills PHS if µ = U/2 and is in general SU(2)-symmetric [63].In the case of zero interaction U = 0, the local generalized charge susceptibility reads where G(iν) = T c † νσ c νσ is the one-particle Green's function.χ νν ′ c is diagonal, hence the eigenvalues can be immediately read from Eq. ( 6).These become doubly degenerate (λ ν = λ −ν = 1/ν 2 ) at PHS, i.e. µ = 0, while they form complex conjugate pairs (λ ν = λ * −ν ) at finite µ.In the left column of Fig. 3 the real (a) and imaginary part (b) of λ i are shown for different µ at finite temperature T = 1/100.
At finite interaction U > 0, χ νν ′ c becomes a more complicated expression [17,[64][65][66], given in the supplemental material [67].The crucial point is the appearance of progressively larger off-diagonal components.The resulting eigenvalues λ i are shown in the right column of Fig. 3. Significantly, at PHS, λ i are still purely real but no longer degenerate.Importantly, this remains true in a finite region of µ around U/2. Far away from PHS, however, the effect of the interaction weakens, and all eigenvalues become complex conjugate pairs.To switch between these two regimes, two eigenvalues have to coalesce: this creates a pair of distinct EPs in µ-space, which delimit and protect the real-eigenvalue "lens"-shaped structure (Fig. 3(c)).In the 2×2-picture of Eq. ( 4), we can identify the interaction U as responsible for the presence of the off-diagonal finite elements c and d in the matrix, and the finite µ for the diagonal element b, which are the two ingredients necessary to satisfy the EP conditions.Hence, for the AL any finite value of U will result in exceptional points away from PHS and a finite-size real eigenvalue lens shape.
Implications on correlation-induced instabilities.We now turn to a more generic scenario, namely the single-orbital Hubbard Hamiltonian on a lattice: with constant hopping t between neighboring sites i and j.This model is again SU(2)-symmetric and for µ = U/2 it fulfills PHS.Except for one or infinite spatial dimensions, the model has not been solved analytically.In order to get a non-perturbative, albeit approximated many-body solution, we use dynamical mean-field theory (DMFT) (which becomes exact only in the limit of infinite dimensions) [68,69] with a continuous-time quantum Monte Carlo solver from w2dynamics [70].As shown in Ref. [30], the eigenvalues λ i and the corresponding eigenvectors v ν i of the local generalized susceptibility ) play an important role for the response functions of the whole lattice: they lead to an enhancement and in some cases to a divergence of the uniform (i.e. for zero transfer momentum q = 0) susceptibility.In particular, for the Bethe lattice with infinite connectivity (where DMFT is exact) the static uniform charge response, obtained by summing the generalized susceptibility over all Matsubara frequencies χ q=0 = νν ′ χ νν ′ q=0 , can be re-expressed in terms of λ i and corresponding weights ).This leads to the following expression which diverges -thus inducing a phase instability in the charge sector -when one eigenvalue fulfills the condition λ i = −1/t 2 .Close to this condition, λ i gives the dominating contribution to the charge response and determines the stability of the physical solution [71].Importantly, this is possible only when λ i is real [72].Although Eq. ( 8) is only exact in the case of the Bethe lattice, numerical calculations have shown [30] that it holds also for a square lattice if t is replaced by a temperature-and µ-dependent t eff (µ, T ).Here, the central role of EPs becomes apparent: their presence guarantees that the imaginary part of λ i remains zero in the whole extended region of the lens shape.In other words, the possibility of inducing a divergence in χ q=0 is not accidental and does not rely on a fine-tuning of U , T and µ: the phase instability is in fact topologically protected.For the square lattice, this is illustrated in Fig. 4, where we plot the real (a) and imaginary part (b) of the eigenvalues λ i of the local charge susceptibility χ νν ′ c close to the critical point of the phase separation (cf.sketch in Fig. 1).Here, the lowest eigenvalue λ I satisfies (up to numerical accuracy) the condition λ I = −1/t 2 eff in the region of the lens shape.Hence, the phase instability condition is also fulfilled for any further reduction of the temperature T (or for any moderate reduction of the interaction U ).In particular, at lower T , we enter a regime, where a first order phase separation occurs.This regime is, thus, characterized by two locally stable DMFT solutions (i.e., two coexisting values of λ I ), corresponding to a less correlated metallic and a "bad metal" phase (connected by an instable solution, where λ I < −1/t 2 eff [71]).Here, the topological robust arguments related to the condition λ I = −1/t 2 eff remain nonetheless applicable, albeit to the two corresponding metastable solutions [73].
Finally, let us notice that a negative eigenvalue is a necessary condition for the instability criterion to be fulfilled [30].Remarkably, the role of the negative eigenvalues in the generalized charge susceptibility has been recently related to the local moment formation [74][75][76][77] and, on a more formal level, to divergences of the irreducible vertex function and the multivaluedness of the Luttinger-Ward functional [17,30,34,56,.Therefore, these negative eigenvalues can be regarded as a feature of strong electronic correlations, which cannot be commonly described by "perturbative" theories, e.g., the random phase approximation.However, the considerations behind Eq. ( 8) are not solely restricted to negative eigenvalues.They can be also applied to the opposite case, where a positive eigenvalue reaching a maximum (e.g., λ i = 1/t 2 eff ) triggers a phase instability, such as the antiferromagnetic transitions of the Hubbard model [99].Thus, in these strongly correlated systems, the presence of EPs is found to generally promote phase instabilities in the ph-channel [100] to a stable phenomenon, and thereby enables the instabilities to naturally occur for a finite range in parameter space.

Conclusion.
We have found the opening of an EP phase for the associated eigenvalues of the static local susceptibility in the U/µ phase diagram of models for correlated electron systems.The remarkable consequence is that the interaction-induced charge instabilities such as the phase separation occurring close to the Mott metal-toinsulator transition in the Hubbard model do not need any fine-tuning but can occur in an entire finite range of parameters.This unexpected global robustness is a consequence of the peculiar centroHermitian form of the susceptibility matrix, which is not dictated by some adhoc antiunitary symmetry but by the time-independence of H and the intrinsic nature of Fermi-Dirac statistics.
The susceptibility EPs represent a clear-cut and compelling manifestation of non-Hermitian topology, surpassing the conventional realizations based on spectral functions.This phenomenon is indeed ubiquitous even in the simplest correlated fermion models and does not require any assumption on the interaction nor any specific choice of non-Hermitian Hamiltonian terms.Our results call for future investigations beyond the local correlation effects on the charge sector considered here: e.g., of the spin or particle-particle channel and including non-local correlations in the description.Further, one could also search for higher order exceptional degeneracies in the susceptibility spectrum and explore the respective consequences on the phase instabilities.This may open new doors to experimentally detectable hallmarks of non-Hermitian topology.

DATA AVAILABILITY
A data set containing all numerical data and plot scripts used to generate the figures of this publication is publicly available on the TU Wien Research Data repository [1].

PROPERTIES OF THE GENERALIZED SUSCEPTIBILITY
The generalized susceptibility in the particle-hole (ph) channel of the main text explicitly reads where c αi e −Hτi .Taking the complex conjugate of this quantity, we get By changing the integration variables from where e iν (′) β = −1, and interchanging the integration limits we obtain the result of the main text.From the definition of the imaginary time-ordering, it immediately follows

LOCAL SUSCEPTIBILITIES WITH SPIN CONSERVATION
In this section, we analyze the matrix properties for the local generalized susceptibilities χ νν ′ ph,σ1...σ4 of a oneorbital model with conservation of spin σ i .This allows us to restrict to the longitudinal (χ νν ′ ph,σσ ′ := χ νν ′ ph,σσσ ′ σ ′ ) and transverse (χ νν ′ ph,σσ ′ := χ νν ′ ph,σσ ′ σ ′ σ ) spin components, that satisfy the following two relations: and are therefore centro-and perHermitian matrices, respectively.PerHermitian matrices [6] are invariant under a transformation that combines complex conjugation with mirror symmetry with respect to the antidiagonal.They belong to the class of κ-Hermitian matrices , where Π refers to any permutation matrix, and have either real or complex conjugate eigenvalues.

CentroHermitian
As stated in the main text, a centroHermitian matrix C has the following property [8,9] with J νν ′ = δ ν(−ν ′ ) .Eigenvalues λ are either real λ = λ * or in complex conjugate (c.c.) pairs ∃λ → ∃λ * .For the eigenvectors Cv = λv, the following relation ∃v → ∃J v * can be easily shown.For individual real eigenvalues, v and J v * are linear dependent From this it follows that also the corresponding weights w of the eigenvalues λ are either real or appear in c.c. pairs: λ ∈ R : λ, λ * ∈ C : PerHermitian A perHermitian matrix P has the following property [6] (P νν ′ ) † = P −ν−ν ′ or equivalently Eigenvalues λ are also either real λ = λ * or in c.c. pairs ∃λ → ∃λ * .Similar for the centroHermitian case for the eigenvectors P v = λv we find the relation Hence, for every right eigenvector v there exists a left eigenvector v † J .For individual real eigenvalues, this left eigenvector is proportional to the inverse of the right eigenvector v −1 = αv † J .From this it again follows that also the corresponding weights w of the eigenvalues λ are either real or appear in c.c. pairs: λ, λ * ∈ C : If a centro-or perHermitian matrix M is further symmetric M = M T [10] we can find an orthonormal eigenbasis v T v = 1 for the subspace where M is diagonalizable.For the non-diagonalizable subspaces v are quasinull v T v = 0.The non-diagonalizable subspace consists of the EPs.(λ ∈ R): For the orthonormal eigenbasis we can show v = x + iy (20) thus α 2 = 1 =⇒ α = ±1.This means that either the real part v is symmetric J x = x and the imaginary part is antisymmetric J y = −y or vice-versa J x = −x and J y = y.This implies a positive weight w > 0 for v = v s with a symmetric real part and negative weight w < 0 for v = v a with an anti-symmetric real part.C.c. pairs λ, λ * ∈ C: For c.c. pairs v ⊥ J v * .Hence, v ̸ = ±J v * and x T J x + y T J y = 0.This is fulfilled for v At an exceptional point (EP) the eigenvector is quasi null v T v = 0 and the eigenvalue λ ∈ R, hence x T x = y T y and v = e iφ J v * .

GENERALIZED CHARGE SUSCEPTIBILITY OF THE ATOMIC LIMIT
Following Refs.[4,11,12] we define the generalized charge susceptibility of the AL in the notation of the main text (with zero bosonic transfer momentum, zero magnetic field, and multiplied with an additional factor of 1/β) as where Z = 1 + 2e βµ + e β(2µ−U ) and n = (e βµ + e β(2µ−U ) )/Z.Note that the term in the second line of Eq. ( 23) becomes proportional to δ νν ′ for PHS at µ = U/2.

pp-CHANNEL
For the local generalized susceptibility χ νν ′ pp,σ1...σ4 of a one-orbital model with conservation of spin σ i the longitudinal particle-particle (pp) channel has quite different spectral properties compared to the ph-channel.χ νν ′ pp,σ1...σ4 is defined as Again taking its complex conjugate and considering spin conservation we obtain two relations: The transversal spin components χ νν ′ pp,σσ ′ in the ppchannel form a perHermitian matrix, however (25a) shows that the longitudinal components χ νν ′ pp,σσ ′ now make up a Hermitian matrix and EPs are not required to ensure real eigenvalues λ for an instability condition.

INSTABILITY CONDITION IN THE REGION OF PHASE SEPARATION
In Fig. 1, we schematically illustrate the fulfillment of the instability condition in the region of the phase sepa-ration (turquoise background).Here, two locally stable DMFT solutions (i.e., two coexisting values of λ I ), corresponding to a less correlated metallic and a "bad metal" phase are connected by an unstable solution (red dotted line), where λ I < −1/t 2 [13].The instability condition λ I = −1/t 2 is fulfilled in the region of the "lens shape" (gray background) between the exceptional points, where λ I is real.For the two locally stable DMFT solutions of λ I , the one which meets the instability criterion corresponds to the metastable solution (the other is thermodynamically stable).
Schematic illustration of the fulfillment of the instability condition in the region of the phase separation (turquoise background).The real parts of the lowest eigenvalue λI and the associated eigenvalue λII , for which λI and λII coalesce at the exceptional points (red dots), are displayed.λI is shown as dotted line below the limit −1/t 2 to indicate the instability of this solution.In the region of the "lens shape" (gray background) both eigenvalues have a zero imaginary part.

GENERALIZED SUSCEPTIBILITY MATRIX PROPERTIES FOR BOSONS
By taking bosonic b ( †) instead of fermionic c ( †) operators for the generalized susceptibility we get where b αi e −Hτ are now the Fourier transforms of bosonic (creation) annihilation operators.The corresponding generalized susceptibility is then describing the propagation (and scattering) of two Bosons instead of two Fermions.The crucial difference of Eq. ( 26) to Eq. (1) in the main text are the bosonic (even) Matsubara frequencies iω (′) = 2n (′) π/β, n (′) ∈ Z.For χ ωω ′ ph,α1...α4 we can also follow the steps in section "Properties of the Generalized Susceptibility" of this Supplemental Material (now e iω (′) β = 1 ) and arrive at the same relation we have found for Fermions: Hence, for Bosons the corresponding χ ββ ′ ph matrix also satisfies the relation where Π ββ ′ is the permutation matrix of β := (ω, α 1 , α 2 ) → β ′ := (−ω, α 2 , α 1 ), and χ ββ ′ ph has the same properties for its eigenvalues λ: they are either real or complex conjugate pairs.However, although the general matrix property remains the same, the choice of fermionic or bososnic Matsubara frequencies does have an important implication on the eigenspectrum of the matrix.The crucial difference is the presence of the zero frequency iω = 0 among the bosonic Matsubara frequencies which is absent for fermionic ones.In addition to the expected complex conjugate eigenvalue pairs and EPs of the fermionic system, the bosonic system has, thus, an extra eigenvalue, which remains always real, even without further symmetries.
Let us then assume for a moment that a similar mechanism, as the one we investigated for the fermionic system, might be also responsible for a phase instability of the bosonic system.Then, since such an instability requires a real eigenvalue to reach a certain threshold, EPs might or might not play a crucial role, however they are no longer a necessary condition for the fulfillment of the instability, due to the extra real eigenvalue.
Consider the simple example of the non-interacting susceptibility of the atomic limit for Fermions where 1 β νν ′ χ νν ′ F = ∂ µ n F (−µ), and for Bosons where /(e βϵ + 1) and n B (ϵ) = 1/(e βϵ − 1)).The resulting eigenvalue spectrum λ i is displayed in Fig. 2. By comparing the spectra of the bosonic (right column) with the fermionc system (left column) an additional eigenvalue λ 0 = −1/µ 2 can be found, which then remains real also out of particlehole symmetry independent of the value of the chemical potential µ.Interestingly in this simple example, λ 0 is the eigenvalue responsible for the (negative) divergence of ∂ µ n B at µ = 0, generically associated to the possible onset of Bose-Einstein condensations.

FIG. 3 .
FIG. 3. Real part (top row) and imaginary part (bottom row) of the eigenvalues λi of χ νν ′ c for the atomic limit (AL) at temperature T = 1/100, U = 0 (a, b) and at U = 1 (c, d) as function of chemical potential away from particle-hole symmetry (PHS) at µ = U/2.The ten eigenvalues which are lowest in Re λi at U = 1 are displayed.

FIG. 4 .
FIG. 4.Real part (a) and imaginary part (b) of the eigenvalues λi of χ νν ′ c for the square lattice Hubbard model (with half bandwidth D = 4t = 1) solved within DMFT as function of chemical potential away from PHS at µ = U/2.Interaction strength U = 2.4 and temperature T = 1/53 coincide with Ref.[30] to show the situation close to the thermodynamic instability.Calculated data are displayed as dots, the positive (µ − U/2)-axis is mapped from the negative one, exploiting the symmetry of the model considered.The ten eigenvalues which are lowest in Re λi at µ − U/2 = 0 are displayed.

FIG. 2 .
FIG. 2. Real part (top row) and imaginary part (bottom row) of the eigenvalues λi of the generalized charge susceptibility for the non-interacting (U = 0) atomic limit at temperature T = 1/100 for Fermions χ νν ′ c (left column) and Bosons χ ωω ′ c (right column) as function of chemical potential.The four resp.five eigenvalues stemming from the inner Matsubara frequencies are displayed.The positive (non physical) values of the chemical potential µ for the bosonic system are included for completeness.