Yukawa-Lorentz Symmetry in Non-Hermitian Dirac Materials

Lorentz spacetime symmetry represents a unifying feature of the fundamental forces, typically manifest at sufficiently high energies, while in quantum materials it emerges in the deep low-energy regime. However, its fate in quantum materials coupled to an environment thus far remained unexplored. We here introduce a general framework of constructing symmetry-protected Lorentz invariant non-Hermitian (NH) Dirac semimetals (DSMs), realized by invoking masslike anti-Hermitian Dirac operators to its Hermitian counterpart. Such NH DSMs feature purely real or imaginary isotropic linear band dispersion, yielding a vanishing density of states. Dynamic mass orderings in NH DSMs thus take place for strong Hubbardlike local interactions through a quantum phase transition, hosting a non-Fermi liquid, beyond which the system becomes an insulator. We show that depending on the internal Clifford algebra between the NH Dirac operator and candidate mass order-parameter, the resulting quantum-critical fluid either remains coupled with the environment or recovers full Hermiticity by decoupling from the bath, while always enjoying an emergent Yukawa-Lorentz symmetry in terms of a unique terminal velocity. We showcase the competition between such mass orderings, their hallmarks on quasiparticle spectra in the ordered phases, and the relevance of our findings for correlated designer NH Dirac materials.

Lorentz spacetime symmetry represents a unifying feature of the fundamental forces, typically manifest at sufficiently high energies, while in quantum materials it emerges in the deep low-energy regime.However, its fate in quantum materials coupled to an environment thus far remained unexplored.We here introduce a general framework of constructing symmetryprotected Lorentz invariant non-Hermitian (NH) Dirac semimetals (DSMs), realized by invoking masslike anti-Hermitian Dirac operators to its Hermitian counterpart.Such NH DSMs feature purely real or imaginary isotropic linear band dispersion, yielding a vanishing density of states.Dynamic mass orderings in NH DSMs thus take place for strong Hubbardlike local interactions through a quantum phase transition, hosting a non-Fermi liquid, beyond which the system becomes an insulator.We show that depending on the internal Clifford algebra between the NH Dirac operator and candidate mass orderparameter, the resulting quantum-critical fluid either remains coupled with the environment or recovers full Hermiticity by decoupling from the bath, while always enjoying an emergent Yukawa-Lorentz symmetry in terms of a unique terminal velocity.We showcase the competition between such mass orderings, their hallmarks on quasiparticle spectra in the ordered phases, and the relevance of our findings for correlated designer NH Dirac materials.
From classical and quantum electrodynamics involving photons to quantum chromodynamics encompassing gluons, field theoretic formulations of fundamental forces rest on the unifying bedrock of the Lorentz symmetry [1,2], typically realized at sufficiently high energies [3][4][5][6][7][8].On the other hand, in Dirac crystals, realized in a number of quantum materials [9][10][11], although such a symmetry may not necessarily be present at the microscopic level, it rises as an emergent phenomenon in the deep infrared regime through boson mediated interparticle interactions [12][13][14][15][16][17].The bosonic degrees of freedom can be vectorlike particles such as helicity-1 photons which interact with Dirac fermions and charged Cooper pairs or spinless scalar order-parameter fluctu-ations at the brink of dynamic mass generation for Dirac quasiparticles via spontaneous symmetry breaking.Irrespective of these microscopic scenarios, the emergent Lorentz symmetry always manifests a unique velocity in the medium, generically tagged as the 'speed of light', not necessarily c.Necessity and elegance of this commonly occurring space-time symmetry often allow us to take it for granted.However, its fate when quantum materials interact with the environment thus far remains an untouched territory, which we here set out to explore.
In the Hamiltonian language, system-to-environment couplings can be modeled by non-Hermitian (NH) operators.But, a one-to-one correspondence between them is still missing.Despite this limitation, we showcase a possible emergent Lorentz symmetry in nonspatial symmetry protected NH Dirac semimetals (DSMs), captured by Lorentz invariant NH Dirac operators, possessing purely real or imaginary eigenvalue spectra.When such an NH DSM arrives at the shore of dynamic mass generation triggered by Hubbardlike finite-range Coulomb repulsion, mediating boson-fermion Yukawa interactions, the resulting strongly coupled incoherent non-Fermi liquid always features a unique terminal velocity in the deep infrared regime.Depending on the pattern of the incipient Dirac insulation, the system achieves the Yukawa-Lorentz symmetry either maintaining its coupling with the environment or by decoupling itself from the bath.These outcomes possibly suggest a generic Lorentz symmetry in NH Dirac materials, despite the variety of its interactions with the environment.
In this work, we establish a general framework for constructing symmetry-protected Lorentz-invariant NH DSMs across various dimensions, described by the Dirac Hamiltonianlike operator, Eq. ( 2), which is achieved by introducing a masslike anti-Hermitian Dirac operator alongside its Hermitian counterpart.Consequently, the thermodynamic, transport, and elastic responses of such NH DSMs closely resemble those of conventional Dirac materials, but in terms of an effective Fermi velocity [Eq.(2)].By employing a leading order ε = 3 − d expansion about d = 3 upper critical dimension of the Gross-Neveu-Yukawa quantum-critical theory, we show that the system in the vicinity of a NH DSM-insulator quantum phase transition features an emergent Yukawa-Lorentz symmetry in terms of a unique velocity, as ex-plicitly shown in Eqs. ( 8)- (12).See also Figs. 1 and 2. Finally, we unfold the quantum (multi-)critical phenomena in correlated NH DSMs (Fig. 3), and characterize it through the critical exponents, such as bosonic and fermionic anomalous dimensions [Eq.( 17)], and correlation length exponent [Eq.(20)].

Minimal model
We set out to construct a minimal effective Hamiltonianlike NH Dirac operator (H NH ), describing a collection of linearly dispersing gapless quasiparticles in d spatial dimensions coupled to environment, such that H NH possesses either purely real or purely imaginary eigenvalues.We begin with the standard Dirac Hamiltonian of the form where v H bears the dimension of the Fermi velocity, Hermitian Γ matrices satisfy the anticommuting Clifford algebra Internal structure of the Dirac spinor Ψ k depends on the microscopic details of the system.For a minimal four-component Dirac spinor, the maximal number of mutually anticommuting Hermitian matrices is five, out of which three (two) can be chosen to be purely real (imaginary).Although our construction is applicable in arbitrary dimensions, here we primarily concentrate on d = 2.Then, without any loss of generality, we choose Γ 1 and Γ 2 to be purely imaginary.
Possible mass terms, Ψ † k M Ψ k , producing isotropically gapped ordered ground states, are then represented by the Hermitian matrices (M ) that anticommute with the Dirac Hamiltonian in Eq. ( 1), namely {M, h 0 } = 0 with M 2 = 1.In d = 2, there are four such mass matrices for a four-component Dirac system M ∈ {M 1 , M 2 , M 3 , Γ 12 }, where Γ jk = iΓ j Γ k .While M j s are purely real for j = 1, 2, 3, each of which breaks the SU(2) chiral symmetry of H D , generated by {M 12 , M 23 , M 31 }, with M jk = iM j M k , the purely imaginary mass matrix Γ 12 transforms as a scalar under the chiral rotation.It breaks the timereversal symmetry under K (complex conjugation).
The crucial observation is that the product of a mass matrix and the Hamiltonian h 0 is anti-Hermitian, namely (M h 0 ) † = −M h 0 .Therefore, we define the NH Dirac operator as a minimal extension of H D [Eq.(1)] by including such a masslike anti-Hermitian term, leading to The real parameter v NH , also bearing the dimension of the Fermi velocity, quantifies the strength of the systemto-environment coupling.The spectrum of the NH Dirac operator real (imaginary NH is the effective Fermi velocity of NH Dirac fermions.Hereafter, we consider the case with the real energy eigenvalues, unless otherwise stated. The form of the NH Dirac operator (H NH ) is restricted by four nonspatial discrete unitary and antiunitary symmetries, among which time-reversal, particle-hole and pseudo-Hermitianity [18], with the corresponding representations given in Table I for a specific choice of M = M 1 .As explicitly shown, once the form of H NH is fixed (for a given M ), none of the four constant mass terms is invariant under all of these symmetries, and therefore cannot be added to H NH without breaking at least one of them.Hence, the nonspatial symmetries protect NH gapless Dirac fermions, irrespective of the specific choice of M and the dimensionality of the system (d).See Methods and Supplementary Note 1 for details.

Scaling
The form of H NH is manifestly Lorentz invariant, with the dynamical exponent z = 1, measuring the relative scaling between the energy (E NH ) and momentum (k).This, in turn, implies that the density of states vanishes in a power-law fashion ρ(E) As a consequence, a weak local (Hubbardlike) short-range interaction is irrelevant for d > 1 and the concomitant quantum phase transition (QPT) into an ordered phase takes place through a strongly coupled quantum critical point (QCP) with its critical behavior described by a Gross-Neveu-Yukawa (GNY) theory, about which more in a moment.Nonetheless, a reduced effective Fermi velocity v F < v H is expected to enhance the mass ordering propensity in NH DSMs in comparison to its counterpart in Hermitian Dirac materials, which we show shortly.
To gain a further insight into the nature of such an NH DSM, next we analyze its response to an external electromagnetic field in terms of the linear longitudinal optical conductivity at zero temperature and finite frequency (ω) within the Kubo formalism.The component of the current operator in the l th spatial direction is

Mass terms
Symmetries Nonspatial symmetries of the non-Hermitian (NH) Dirac operator.We here show the symmetries of the NH operator HNH in d = 2 spatial dimensions with the anti-Hermitian term containing the matrix M = M1 [Eq.( 2)].T+ and C+ (T− and C−) represent time-reversal (antiunitary particle-hole) symmetry, while PH1,2 (PSH1,2) are unitary particlehole (pseudo-Hermiticity) symmetry, under which . We work with a representation in which Hermitian matrices M1, M2 and M3 (Γ1 and Γ2) are purely real (imaginary).The symbol ✓ (×) indicates whether a mass term is even (odd) under a specific symmetry.Notice that none of the four possible uniform mass terms is invariant under all the nonspatial symmetries of HNH, and thus can not be added to HNH without breaking at least one of them.A similar symmetry protection for NH Dirac nodes in d = 2 exists for M = M2, M3 and iΓ1Γ2 (purely imaginary), and also in three dimensions.See Methods and Supplementary Note 1 for details.
) is the fermionic Green's function.After the analytic continuation iω → ω + iδ with δ > 0 to real frequency ω, and using the Kubo formula, we obtain the optical conductivity δ lm and σ (3) in d = 2 and d = 3, respectively.Here N f is the number of four-component fermion flavors.See Methods and Supplementary Note 2. Therefore, NH DSMs and its Hermitian counterparts share the same value of the optical conductivity [19][20][21], with the Fermi velocity in d = 3 replaced by the effective one for NH Dirac fermions (v F ). Furthermore, the frequency-dependent optical shear viscosity of noninteracting NH DSMs in d dimensions features a single independent component due to the underlying rotational symmetry, and is given by η Therefore, NH and Hermitian DSMs possess the same value of the optical shear viscosity [22], with the Fermi velocity being replaced by an effective one in the former systems.See Methods and Supplementary Note 3.These outcomes may imply that the associated conformal field theories feature the same central charge.

Mass orders
The order parameter (OP) of a uniform Dirac mass that isotropically gaps out Dirac fermions can be expressed as an n-component vector But, the appearance of a mass matrix (M ) in H NH frag-ments the mass OPs and they can be classified according to the canonical (anti)commutation relations between N j s and M : (a) commuting class mass (CCM) for which [N j , M ] = 0 for all j, (b) anticommuting class mass (ACM) with {N j , M } = 0 for all j, and (c) mixed class mass (MCM) for which Hereafter, we mainly focus on the former two classes of mass ordering.
The ordering tendencies toward the nucleation of these two classes of mass ordering, revealing the impact of the NH parameter (v NH ), can be estimated from their corresponding (nonuniversal) bare mean-field susceptibilities at zero external frequency and momentum, given by respectively, where ), and Λ is the ultraviolet momentum cutoff up to which the energy-momentum relation remain linear.See Supplementary Note 4. As the effective Fermi velocity in an NH DSM (v F ) is smaller than its counterpart in Hermitian Dirac materials (v H ), the bare mean-field susceptibilities are larger in the former system.Given that the (nonuniversal) critical coupling constant (u ⋆ ) for any mass ordering is inversely proportional to the bare mean-field susceptibility, NH DSMs are conducive to mass formation at weaker interactions, stemming from its enhanced density of states.However, the competition between CCMs and ACMs requires notion about their associated quantum critical behaviors, captured from the renormalization group (RG) flows of the velocity parameters, v H , v NH and v F , which we discuss next.This analysis will also shed light on the emergent multi-critical behavior near the MCM condensation.Notice that the dimen- sionless susceptibility (χv F /Λ d−1 ) and critical interaction strength (u ⋆ /[v F Λ d−1 ]) are cut-off independent and universal.In fact, the later one can be compared with the critical interaction to the band width ratio in a microscopic setup, which for the Néel antiferromagnet order in Hermitian honeycomb lattice is U ⋆ /t ≈ 4.5 ± 0.5 [23], where t (U ) is the nearest-neighbor hopping amplitude setting the Fermi velocity (on-site Hubbard repulsion).

Quantum critical theory
The dynamics of the bosonic OP fluctuations in a GNY quantum critical theory, describing a strong-coupling instability of NH DSMs, is captured by the Φ 4 theory, with the bosonic correlator given by G , where v B is the velocity of the OP fluctuations.The dynamics of fermions is governed by the NH Dirac operator in Eq. ( 2).These two degrees of freedom are coupled through a Yukawa vertex, entering the imaginary time (τ ) action as At the upper critical three spatial dimensions, both Yukawa (g) and the Φ 4 coupling (λ) are marginal [24].We, therefore, use the distance from the upper critical dimension ε = 3 − d as the expansion parameter to capture the low-energy phenomena close to the GNY QCP.The ultraviolet divergences of the Feynman diagrams are captured using the dimensional regularization and the cut-off independent universal physical quantities, such as the RG flow equations and critical exponents, are computed using the method of minimal subtraction [2,24].In particular, we are interested in the fate of an emergent Yukawa-Lorentz symmetry close to such a RG fixed point, since fermionic and bosonic velocities are generically different at the lattice (ultraviolet) scale.

Yukawa-Lorentz symmetry
To this end, we compute the fermionic [Fig.1(a)] and bosonic [Fig.1(b)] self-energy diagrams to the leading order in the ε expansion.Irrespective of the nature of the mass ordering (CCM or ACM), the RG flow of the Hermitian component of the Fermi velocity (v H ) takes the form after rescaling the Yukawa coupling according to g 2 → g 2 /(8π 2 ), where β Q ≡ dQ/d ln b and b is the RG time.
The flow equations for the remaining two velocities, when the mass matrix commutes with M (CCM) are and These RG flow equations [Eq.( 8)- (10)] imply that at the GNY QCP with g 2 ∼ ε, the terminal velocities are such that , with all three velocities being nonzero, independent of their initial values.See Supplementary Note 5. We also confirm this outcome by numerically solving these flow equations.See Fig. 2(a).Therefore, a new fixed point with an enlarged symmetry emerges, at which the system remains coupled to the environment and the effective Fermi velocity of NH Dirac fermions (v F ) is equal to the bosonic OP velocity (v B ), with the nonzero terminal values for v H , v NH and v B .We name it non-Hermitian Yukawa-Lorentz symmetry.
On the other hand, when a NH DSM arrives at the brink of the ACM condensation, the flow equations for v NH and v B take the forms (see Supplementary Note 5) and respectively.These two flow equations along with Eq. ( 8) imply that near the ACM class GNY QCP (g 2 ∼ ε), the system recovers Hermiticity by decoupling itself from the environment and a conventional Yukawa-Lorentz symmetry emerges with v NH = 0 and v F = v H = v B , irrespective of their bare values.We further confirm this outcome by numerically solving the flow equations.See Fig. 2(b).

Mean-field theory
Guided by the emergent quantum critical theory, we now compare the tendencies between CCM and ACM orderings.At the corresponding GNY QCPs, their renormalized susceptibilities are respectively  19).The green fixed point in (a) with g 2 1 ̸ = g 2 2 controls a first order transition between the CCM and ACM, which can be tuned by one of the two bosonic masses, and thus devoid of any single-parameter scaling.Their distinct RG flow equations are shown in Eq. (19).By contrast, green dot in (b) with g 2 1 = g 2 2 and an enlarged O(n) symmetry drives a continuous DSM-mixed-class-mass QPT, tunable by a single mass parameter m 2 = m 2 1 = m 2 2 with its RG flow equation shown in Eq. ( 21).obtained from Eq. ( 6) by replacing various velocities with their terminal ones.As is the critical strength of the coupling constant for the CCM (ACM) ordering.This outcome is particularly important when the mass matrix M , appearing in H NH [Eq.( 2)], is a component of a vector mass order of a Hermitian Dirac system.Any such choice of M immediately fragments the vector mass order into CCM and ACM, and our susceptibility analysis suggests that the OP corresponding to the CCM enjoys a stronger propensity toward nucleation.
This prediction can be tested from the quasiparticle spectra of NH Dirac fermions inside mass ordered phases, described by an effective NH single-particle operator where N j s are the Hermitian mass matrices, satisfying {N j , h 0 } = 0 for all j, and ∆ j is the OP amplitude.The excitation spectrum of H MF NH crucially depends on the nature of the mass ordering.Namely, (a) for [N j , M ] = 0 (CCM), the energy eigenvalues are either purely real (v H > v NH ) or complex (v H < v NH ), whereas (b) for {N j , M } = 0 (ACM) the excitation spectrum is complex for any v H and v NH .Next we exemplify these outcomes by focusing on a specific microscopic scenario.

NH honeycomb Hubbard model
Consider a collection of massless Dirac fermions on graphene's honeycomb lattice.In this system, the Neél antiferromagnetic (AFM) order corresponds to the threecomponent vector mass that breaks the O(3) spin rotational symmetry, favored by a strong on-site Hubbard repulsion at half-filling [23] due to the underlying bipartite nature of the lattice allowing the formation of staggered magnetization without any frustration.We construct an NH Dirac operator by choosing the easy z-axis component of the AFM order as M , which then naturally becomes a CCM order.By contrast, two planar or easy-plane components of the AFM order fall within the category of ACM.See Supplementary Note 6.Our susceptibility calculation predicts that such an NH Hubbard model should prefer nucleation of an Ising symmetry breaking easy-axis AFM order over the XY symmetry breaking easy-plane one, and the quasi-particle spectra inside the ordered phase must be purely real unless v H < v NH at the bare level.These predictions can serve as the litmus tests for our quantum critical theory of correlated NH Dirac materials in quantum Monte Carlo simulations and exact numerical diagonalizations [25].

NH criticality
Finally, we compute the critical exponents near the GNY QCPs associated with the n 1 -component CCM and n 2component ACM orderings.The RG flow equations for the corresponding Yukawa couplings respectively read as after rescaling the coupling constants according to See Supplementary Note 7 for details.The Yukawa fixed points at g 2 [blue dot in Fig. 3(a)] respectively control the continuous QPTs into the CCM and ACM orders.At the Yukawa fixed points the fermionic and bosonic anomalous dimensions are respectively for j = 1, 2, and the fermionic Green's functions scale as , indicating the onset of non-Fermi liquids therein.The bosonic Green's function scale as . While the universal critical exponents are independent of the velocity parameters, their RG flows shown in Fig. 2 determine the scaledependent (b) fermionic and bosonic Green's or spectral function within the quantum critical regime.Here b can be identified as b = Λ/Λ 0 , where Λ(Λ 0 ) is the running (fixed lattice or ultraviolet) momentum or energy scale.
The last terms in Eqs. ( 15) and ( 16) are pertinent in the proximity to an n-component MCM ordering, where n = n 1 + n 2 , and they are nontrivial only when the bare value of the NH Fermi velocity (v 0 NH ) is zero, corresponding to a Hermitian Dirac system.Then a fully stable quantum multi-critical point with enlarged O(n) symmetry emerges at g 2 , 27].These terms, however, do not exist in NH Dirac system as CCM and ACM components of the MCM order possess distinct Yukawa-Lorentz symmetries (hence cannot be coupled).The fixed point located at (g ) thus controls a direct first-order phase transition between them, as the transition can be tuned by either one of the bosonic masses m 2 1 or m 2 2 , for which the distinct RG flow equations are shown below in Eq. ( 19).We do not discuss such firstorder transition between two ordered phases in any further detail.These two scenarios are shown in Fig. 3.
The RG flow equations for the Φ 4 couplings associated with the CCM and ACM orderings assume the forms for j = 1 and 2, respectively, in terms of the rescaled cou- The fixed point values λ j,⋆ (somewhat lengthy expression) are obtained from the solution of β λj = 0.The RG flow equation for the bosonic mass parameters that serve as the tuning parameter for the NH DSM to CCM and ACM QPTs respectively take the forms for j = 1 and 2, yielding the correlation length exponents On the other hand, the RG flow equation for the mass parameter controlling the continuous QPT between two ordered phases across the multi-critical point [green dot in Fig. 3(b)] in Hermitian Dirac system reads as where

Conclusions
In this work, we develop a general formalism of constructing symmetry protected Lorentz invariant NH DSMs in any dimension by introducing masslike anti-Hermitian Dirac operator to its Hermitian counterpart, featuring purely real or imaginary eigenvalue spectra.Their thermodynamic, transport and elastic responses closely mimic the ones in conventional Dirac materials, however in terms of an effective Fermi velocity of NH Dirac fermions, where v H (v NH ) is the Fermi velocity associated with the Hermitian (anti-Hermitian) component of the NH Dirac operator.Following Eq. ( 2), our construction of the NH DSMs can be generalized to any semimetal of arbitrary energy dispersion relation, described by an Hermitian operator (h 0 ) and accommodates mass orders (M ), to construct its NH counterpart, since in all these cases the anti-Hermitian operator M h 0 exists.We show that when NH DSMs arrive in the close proximity to any mass ordering via spontaneous symmetry breaking, the emergent non-Fermi liquid features NH (for CCM) or Hermitian (for ACM) Yukawa-Lorentz symmetry in terms of a unique terminal velocity for all the participating degrees of freedom.We also determine the nontrivial critical exponents associated with the NH DSM to Dirac insulator QPTs, revealing their non-Gaussian nature in d = 2, favored by strong Hubbardlike local interactions.Combining mean-field theory and RG analysis, we address the competition among different classes of mass orderings and argue that the nature of the Dirac insulators can be identified from the quasiparticle spectra inside the ordered phases, which can serve as a direct test of our proposed effective description of correlated NH DSMs in numerical simulations [25] and experiments.In three spatial dimensions, the NH DSM to insulator QPTs are Gaussian in nature as they take place at g 2 j = λ j = 0 since ε = 0 therein, yielding mean-field critical exponents ν = 1/2 and η Ψ = η Φ = 0. Still, the fermionic and bosonic quasiparticle poles suffer logarithmic corrections, producing a marginal Fermi liquid, which also manifests the emergent Yukawa-Lorentz symmetry.
Even though the universal critical exponents are independent of any velocity parameters and microscopic details of the system, and only depends on the parameter ε, the nature of the order phases (characterized by the number of bosonic order-parameter components n), and the four-component Dirac fermionic flavor number N f , the signatures of v H , v NH and v B can be probed from the energy dependent measurement of the two-point fermionic (G F ) and bosonic (G B ) correlation functions, capturing the RG flow of these parameters shown in Fig. 2. Furthermore, the non-Hermiticity can also change the nature of the ordered phases in NH DSMs in comparison to its counterpart in Hermitian DSMs and that way directly affecting the associated quantum critical phenomena.For example, Hubbard repulsion in Hermitian honeycomb lattice supports O(3) symmetry breaking Néel antiferromagnet order, with the critical exponents determined by setting N f = 2 and n = 3.If we choose the mass matrix M to be the z or easy-axis component of the Néel order in the construction of the NH Dirac operator H NH [Eq.( 2)], the ultimate ordered state is expected to be the Ising symmetry breaking easy-axis antiferromagnet.Therefore, the critical exponent for such an NH honeycomb Hubbard model is set by N f = 2 and n = 1.
For experimental realizations of our proposal, consider a paradigmatic Dirac system in d = 2, graphene [9], where conventional Dirac Hamiltonian (h 0 ) results from electronic hopping between the nearest-neighbor sites of A and B sublattices.As such graphene can accommodate a plethora of Dirac masses [28,29] and any one of them can be chosen as M in H NH [Eq.( 2)].Consider the simplest one, charge-density-wave with staggered pattern of electronic density between two sublattices [30].Then the anti-Hermitian operator M h 0 also corresponds to nearest-neighbor hopping on the honeycomb lattice.But, the hopping strength from A to B sites is stronger or weaker than that along the opposite direction, yielding non-Hermiticity.See Supplementary Note 8.Such a simple NH Dirac operator can be engineered in designer electronic and optical honeycomb lattices, on which Hermitian graphene has already been realized [31,32].In light of the recent proposal of realizing NH one-dimensional optical lattice with a left-right hopping imbalance [33], a hopping imbalance between the A → B and B → A directions can be engineered with two copies of the honeycomb lattice, occupied by neutral atoms living in the ground state and first excited state, which are coupled by running waves along three nearest-neighbor bond directions and the sites constituted by the excited state atoms undergoes a rapid loss.When the wavelength of the running wave is equal to the lattice spacing, a NH Dirac operator on optical honeycomb lattice with hopping imbalance between A → B and B → A directions can be realized.We point out that on one-dimensional optical lattices NH models with coupling to the environment has recently been demonstrated [34].A similar construction can possibly be engineered on designer electronic lattices to mimic the desired hopping imbalance.
Given that Hubbard repulsion driven AFM ordering has been observed on optical graphene [32] and the tunable NH coefficient v NH reduces the critical interaction for this ordering, since AFM is a CCM when M is the charge-density-wave order, our predicted mass nucleation at weak coupling and the associated quantum critical phenomena should be observed in these NH designer Dirac materials.Simplicity of this construction should also make the proposed phenomena observable at least on three-dimensional NH optical Dirac lattices, which also supports a number of mass orders [35].
Altogether our Lorentz invariant construction of the NH Dirac operator constitutes an ideal theoretical framework to extend the realm of various exotic many-body phenomena in the presence of system-to-environment interactions.Among them quantum electrodynamics, topological defects in the ordered phases, magnetic catalysis, superconductivity, chiral anomaly, to name a few in NH Dirac materials, are the prominent and fascinating ones.The present work lays the foundation of systematic future investigations of these avenues.

Methods
In this section, we outline the key methodologies employed in this work to arrive at various conclusions.

Symmetry analysis of the NH Dirac operator
The symmetry analysis of the NH Dirac operator, given by Eq. ( 2) is straightforwardly performed, by using the corresponding canonical (anti)commutation relations of the 16 Hermitian 4 × 4 Dirac matrices squaring to one, in both d = 2 and d = 3 spatial dimensions.For this analysis, it is also of crucial importance to recall that out of these 16 Dirac matrices maximally five of them mutually anticommute.In this set of five mutually anticommuting matrices, three are purely real while two are purely imaginary.Now, the structure of the purely Hermitian Dirac Hamiltonian depends on the dimensionality.In d = 2, the two Dirac matrices in the kinetic term can be chosen to be purely imaginary, while in d = 3 time-reversal symmetry enforces the three Dirac matrices in the kinetic to be purely real.In either case, one of the remaining Dirac matrices can be chosen to form the anti-Hermitian piece of the NH Dirac operator.Given the above, we classify the NH Dirac Hamiltonian according to the possible symmetries of the NH operators [18].This procedure leads to the classification presented in Table I for NH Dirac operator in d = 2, while the analogous classification for the d = 3 is given in Supplementary Note 1, together with the additional details, summarized in Supplementary Table 1 and Supplementary Table 2.

Zero temperature optical conductivity
The optical conductivity at zero temperature is computed using the standard Kubo formula, which, for the completely isotropic system, considered here, takes the form, where Π(iω) is the diagonal element of the polarization tensor of the noninteracting NH Dirac fermions after subtracting its zero-frequency part, Π 1(a).Notice that σ(ω) is computed at real frequency ω after performing the analytical continuation from the Matsubara frequency iω → ω + iδ in the limit δ → 0.
The polarization tensor explicitly reads as where j = 1, • • • , d, the current operator J j = (v H + v NH M )Γ j , and is the fermionic Green's function corresponding to the NH Dirac operator in Eq. ( 2) with v 2 F = v 2 H − v 2 NH .The rotational symmetry implies isotropy of the polarization tensor.The straightforward calculations then yield the results in Eq. ( 4) for the zero-temperature optical conductivity in d = 2 and d = 3.Additional details are given in Supplementary Note 2.

Zero temperature shear viscosity
In the case of the optical shear viscosity tensor, we use the corresponding Kubo formula, given by in terms of the stress tensor correlation function Here, T lm (k) = T

RG flow equations in the GNY theory
We first outline the details of the RG analysis leading to the form of the RG flow equations of the fermionic and bosonic velocities, given by Eqs. ( 8)- (12).To this end, we employ the leading-order ε expansion about d = 3 upper critical spatial dimensions of the NH Gross-Neveu-Yukawa (GNY) theory with ε = 3 − d.We compute the self-energy diagrams for the fermions and bosons, shown in 1(a) and (b), respectively, while the Yukawa bosonfermion interaction is given by Eq. (7).See also Supplementary Figures 1(b) and (c).The fermionic self-energy takes the form while the bosonic self-energy is The fermionic Green's function is given by Eq. ( 24), while the bosonic propagator is We then employ the minimal subtraction scheme [24], as detailed in Supplementary Note 5, to find the flow equations of the fermionic and bosonic velocity parameters, both when the order parameter is CCM and ACM, ultimately yileding Eqs. ( 8)- (12).Additional technical details are given in Supplemental Note 5.
Analogously, we find the flow equations of the Yukawa and the Φ 4 -vertex, by computing the corresponding diagrams, respectively, shown in Supplementary Figure 2(a) and Supplementary Figures 2(b) and (c).We then employ the minimal subtraction scheme to extract the corresponding renormalization constants, which ultimately yield the RG flow equations in Eqs. ( 15), (16), and (18).Additional details are shown in Supplemental Note 6.

FIG. 3 . 2 1 and g 2 2
FIG. 3. Schematic renormalization-group (RG) flow in the Yukawa plane.The flow for (a) non-Hermitian and (b) Hermitian Dirac systems [Eqs.(15) and (16)].The Yukawa couplings g 2 1 and g 2 2 are always measured in units of the RG parameter ε.Colored dots correspond to RG fixed points.The [blue] fixed points in both (a) and (b) control continuous quantum phase transition (QPT) from Dirac semimetal (DSM) to commuting class mass (CCM) [anticommuting class mass (ACM)] order, tuned by the bosonic mass m 2 1 [m 2 2 ].The flow equations of the bosonic masses are shown in Eq. (19).The green fixed point in (a) with g 2 1 ̸ = g 2 2 controls a first order transition between the CCM and ACM, which can be tuned by one of the two bosonic masses, and thus devoid of any single-parameter scaling.Their distinct RG flow equations are shown in Eq.(19).By contrast, green dot in (b) with g 2 1 = g 2 2 and an enlarged O(n) symmetry drives a continuous DSM-mixed-class-mass QPT, tunable by a single mass parameter m 2 = m 2 1 = m 2 2 with its RG flow equation shown in Eq. (21).
lm are the orbital and the spin parts of the stress tensor, respectively, given byT (o) lm (k) = k m ∂H NH (k) ∂k l , T(s)lm (k) = i[S lm , H(k)],(27)with S lm = (i/8)[Γ l , Γ m ] representing the generators of the spin rotations, and i, j, k and l are the spatial indices.In particular, in d = 2, the generator of spin rotations takes the form S 12 = (i/4)Γ 1 Γ 2 , while in d = 3, the three generators can be explicitly written as S lm = (i/4)Γ l Γ m , with l ̸ = m, and l, m = 1, 2, 3.The straightforward calculations then yield Eq. (5).The additional details are given in Supplemental Note 3.