Emergent symmetries and coexisting orders in Dirac fermion systems

Emilio Torres, Lukas Weber, Lukas Janssen, Stefan Wessel, and Michael M. Scherer Institut für Theoretische Physik, Universität zu Köln, 50937 Cologne, Germany Institut für Theoretische Festkörperphysik, JARA-FIT and JARA-HPC, RWTH Aachen University, 52056 Aachen, Germany Institut für Theoretische Physik and Würzburg-Dresden Cluster of Excellence ct.qmat, Technische Universität Dresden, 01062 Dresden, Germany

The quantum phase diagram and critical behavior of two-dimensional Dirac fermions coupled to two compatible order-parameter fields with O(N1) ⊕ O(N2) symmetry is investigated. Recent numerical studies of such systems have reported evidence for non-Landau-Ginzburg-Wilson transitions and emergent O(N1 + N2) symmetry between the two ordered states, which has been interpreted within a scenario of deconfined quantum criticality in (2+1)-dimensional Dirac materials. Here, we provide two theoretical approaches to refine the phase diagrams of such systems. In the immediate vicinity of the multicritical point between the ordered phases and the semimetallic phase, we employ a non-perturbative field-theoretical analysis based on the functional renormalization group. For the particular case of N1 = 3, N2 = 1, we perform a large-scale quantum Monte Carlo analysis of the strong-coupling region, where both orders meet. Our findings support the robust emergence of enhanced symmetry at the multicritical point and suggest the transition between the two ordered phases to take place via a sequence of continuous transitions. In particular, we find that intermediate regimes of coexistence are present in the phase diagram for all values of N1 and N2.
Indeed, recent quantum Monte Carlo (QMC) simulations of Dirac fermion systems [19,30,31] suggest continuous non-LGW transitions between two ordered phases, reminiscent of DQCPs. In particular, the findings in Ref. 30 for a fermionic model on the honeycomb lattice indicate that a system of Dirac fermions with anticommuting masses that break an O(3) and Z 2 symmetry, respectively, supports a line of continuous transitions that separates the two phases, featuring an emergent O(4) symmetry. In particular, no definite signs of coexisting orders were reported in Ref. 30.
Here, we examine the case of a general system of Dirac fermions coupled to two compatible order parameters (OPs) with O(N 1 ) ⊕ O(N 2 ) symmetry by following two different and complementary routes: a non-perturbative field-theoretical renormalization group (RG) approach, i.e., the functional RG (FRG) [32], and a refined QMC analysis for the model in Ref. 30. The non-perturbative FRG can be performed directly in 2 + 1 dimensions and allows us to assess the multicritical behavior of the model more precisely than leading-order expansions [33][34][35]. We firmly establish the emergence of O(N 1 + N 2 ) symmetry at the multicritical point for all consistent values of N 1 and N 2 . Furthermore, our approach facilitates a study of the phases with broken symmetry and we find robust indications for an intermediate coexistence phase for all choices of N 1 , N 2 . Further evidence of coexistence for N 1 = 3, N 2 = 1 is provided by our refined QMC analysis.
The O(N 1 ) ⊕ O(N 2 ) symmetric system defined in Eq. (1) includes a symmetry-enlarged O(N 1 + N 2 )invariant subspace, when u 1 = u 2 = u 3 , m 2 1 = m 2 2 , and g 2 1 = g 2 2 . The leading-order -expansion analysis [33,34] studied the RG fixed point with the enlarged O(N ) symmetry, referred to as the isotropic fixed point (IFP). Within that approach, the IFP is found to be stable for all consistent values of N . However, it is well known from studies of purely bosonic O(N 1 ) ⊕ O(N 2 ) models [8] that the leading-order -expansion severely overestimates the stability of the IFP. In contrast, the FRG provides more faithful results already at the low truncation orders that we exploit here [11].
An important subtlety in Dirac fermion systems concerns the determination of the nature of a multicritical point. In a purely bosonic theory, a multicritical point can either be bicritical or tetracritical, and the two cases can be distinguished by the sign of the quantity ∆ = u 1 u 2 − u 2 3 in terms of the quartic couplings u i at the RG fixed point, which is bicritical if ∆ < 0 and tetracritical if ∆ > 0 [2,8,11]. For the symmetry-enhanced case, ∆ = 0, the above classification is valid if the submanifold in coupling space determined by ∆ = 0 is closed under the RG flow. This is not generally the case in theories with massless fermions as can be shown using the -expansion [33]: the submanifold in theory space defined by u 1 = u 2 = u 3 =: u = 0 satisfies ∆ = 0, but its dependence on the logarithmic RG scale t is given by ∂ t ∆ = 2u(g 2 1 − g 2 2 ) 2 . Therefore, in the presence of Dirac fermions, the sign of ∆ may be subject to change.
Functional renormalization group. We employ the non-perturbative FRG approach [32,44] to evaluate the generating functional of one-particle irreducible n-point correlation functions Γ (n) . See Refs. [45][46][47][48][49][50][51][52][53][54][55] for applications to low-dimensional GNY systems. This method allows us to calculate directly in D = 2 + 1 integrate out the flow also within the ordered phases. Central to the FRG method is the exact renormalization group flow equation is the second functional derivative with respect to all field degrees of freedom, and R k is an infrared cutoff function. This flow equation interpolates between the bare action at the ultraviolet (UV) cutoff k = Λ and the full quantum effective action Γ = Γ k=0 .
We employ a truncation based on the original form of the microscopic action in Eq. (1), i.e., which is known as the extended local potential approxi- mation, LPA . The scale dependence of the Yukawa couplings g 1,k and g 2,k as well as the wavefunction renormalizations Z Φ,k have been made explicit. The bosonic potential V k depends on the O(N 1 ) and O(N 2 ) invariant quantities ρ φ = φ a φ a /2, ρ χ = χ b χ b /2, and we expand its dimensionless form u = k −D V k as a function of the dimensionless O(N i )-symmetric field variables ρ φ = Z φ,k k 2−D ρ φ andρ χ = Z χ,k k 2−D ρ χ around a (running) minimum, denoted κ φ/χ =ρ φ/χ,min . A truncation up to order n max in the fields will be referred to in the following as LPA n max . The nature of the minimum corresponds to three different scenarios, depending on whether the individual OP fields are in the symmetric or in the spontaneously symmetry broken regime: (1) Both OP fields remain in their symmetric phases, i.e., the minimum of the potential is at κ φ/χ = 0. (2) The OP φ acquires an expectation value, while the other OP χ remains in its symmetric phase. Similarly, for the reversed situation.
(3) The minimum of the potential is at nonzero expectation values for both OPs, i.e., κ φ/χ = 0. The FRG flow equations for the dimensionless expansion coefficients of the potential u, including the quadratic mass terms and the quartic couplings, the dimensionless Yukawa couplingsg 2 i,k = g 2 i k D−4 Z −1 i,k Z −2 ψ,k , as well as the anomalous dimensions ∂ t Z ψ/φ/χ are obtained by taking the corresponding derivatives of the ansatz in Eq. (2). Flow equations for the expectation values κ φ/χ are obtained from the condition that they are a minimum, cf. the Supplemental Material [56] for details. Essentially, this procedure yields a set of coupled flow equations of the form Stability of the symmetry-enhanced fixed point. A fixed point is defined by the condition ∀i : Here, we focus on the IFP, facilitating the symmetry- directions, attracted to the fixed point. The IFP is stable, if the third-largest eigenvalue of the stability matrix is negative (called the stability exponent).
We report the numerical results for the five largest eigenvalues of the stability matrix in Tab. I, exhibiting the stability of the IFP for all consistent choices of N . In particular, the stability exponent has a sizable magnitude of O(1). Therefore, any perturbation of the enhanced symmetry near the multicritical point will die out rather quickly, supporting a strong tendency towards the emergent symmetry. This strong tendency towards emergent O(N 1 + N 2 ) is consistent with the QMC findings in Ref. 30. Noticeable, the symmetry-enhancement is naturally realized in our extended LGW approach through the fluctuations of massless Dirac fermions, i.e., without requiring the inclusion of additional topological terms. This contrasts to the case of purely bosonic O(N 1 ) ⊕ O(N 2 ) models [11,57], for which symmetryenhancement is supported only for N 1 = N 2 = 1, with an almost marginal stability exponent.
Phase diagram from field theory. The phase diagram finally can be obtained from integrating the FRG flow equations towards the infrared. To this end, we formulate the initial value problem at an arbitrary ultraviolet scale Λ, used to set the units. To resolve the phase diagram in the vicinity of the symmetry-enhanced quantum multicritical point, perform a sweep of initial conditions in the vicinity if IFP. We consider three different choices for N 1 and N 2 , (i) two coupled Ising order parameters . The resulting phase diagrams are shown in Fig. 1, and we find that all cases exhibit extended coexistence regions, at which both OPs develop a vacuum expectation value.
The lack of data points in the region close to the line δm 2 1 = δm 2 2 in Fig. 1 emerging from the IFP is related to a numerical instability of the FRG flow. In this region, the expectation values κ i,k do not converge to a definite value (either zero or nonzero) in the limit k → 0. This behavior occurs along with the appearance of Goldstone modes in each of the adjacent phases, whose interplay with the massive modes has the tendency to drive the system out of the symmetry-broken phase. Indeed, when the adjacent symmetry-broken phases involve no massless modes (for N 1 = N 2 = 1), the line of exact O(N ) symmetry can be resolved effortlessly, while the region of numerical instability grows with the number of available Goldstone modes. The missing regions of the phase diagrams could be determined by FRG methods beyond the scope of this paper, e.g., by applying pseudo-spectral methods [51,52]. However, the LPA truncation already clearly establishes the appearance of a coexistence region near the IFP.
Quantum Monte Carlo. For the case of O(3) ⊕ Z 2 , we obtain direct support for the coexistence region also from a refined QMC analysis for the microscopic model of Ref. 30. For this purpose, we derive an effective quantum spin model that emerges in the strongly-interacting regime of the interacting Dirac fermion model, and which can be simulated by more efficient QMC methods using cluster updates and larger lattices than accessible to the QMC approach for the original fermionic model. The Hamiltonian of this effective quantum spin model, obtained by perturbation theory about the strong-interaction limit [58], reads in terms of Heisenberg S = 1/2 spins S i , residing on a honeycomb lattice and coupled via bond-centered strengths χ ij to a transverse-field (h) Ising model of spins σ ij , located on the nearest-neighbor bonds of the honeycomb lattice. The summation over i, j ( ij, kl ) extends over nearest-neighbor Heisenberg (Ising model) spins. In terms of the nearest-neighbor hopping t, the Hubbard interaction U , and the fermion-spin couplings ξ ij of the underlying fermionic model, we obtain J ij = 4(t 2 + ξ 2 ij )/U , and χ ij = 8tξ ij /U in second-order perturbation theory [56]. These relations can also be used to specify parameter values of H and compare to the results of Ref. 30.
With a staggered pattern of ξ ij = +ξ, 0, −ξ (cf. the inset of Fig. 2) as in Ref. 30, the model has a combined lattice-inversion and Ising model spin-flip Z 2 symmetry in addition to the O(3) [SU (2)] symmetry of the Heisenberg exchange. To connect to the previous results, we also fix ξ = 0.5, J I = 1, and t = 1, and probe the strongcoupling regime for values of U > 6, where large singleparticle gaps prevailed. We used a hybrid QMC parallel tempering scheme [56,59,60] for the Hamiltonian H on periodic lattices with N H Heisenberg spins (and N I = 3N H /2 Ising spins), for N H up to 2400, based on the stochastic series expansion approach [61][62][63]. In particular, we monitored the evolution of the ferromagnetic Ising OP m I = | 1 NI NI b=1 σ z b | and the antiferromagnetic  h. As an example, Fig. 3 shows the finite-size scaling of both OPs, for h in the transition region at U = 7. The algebraic behaviours at the order-to-disorder transitions are in accord with the anticipated asymptotic scalings for the purely bosonic Ising and Heisenberg universality classes in dimension D = 2 + 1, and yield two distinct critical field strengths, h I c = 3.2088(5) for m I , and h H c = 3.202(2) for m H . The resulting phase diagram in the above parameter regime is shown in Fig. 2. It exhibits a gosshamer intermediate phase of coexisting orders between the small (large) h phase with pure Z 2 (O(3)) symmetry breaking: For small h, the Ising model spins order ferromagnet-ically, spontaneously breaking the Z 2 symmetry. Due to the staggered ξ ij -coupling, this induces a preferred dimerization pattern on the honeycomb lattice, which leads to a dominant valence-bond singlet formation of the Heisenberg spins along the stronger (ξ ij σ z ij > 0) bonds. Upon increasing h, the ferromagnetic order reduces, in effect weakening also the induced dimerization, so that antiferromagnetic Heisenberg spin order can eventually set in prior to the full suppression of the ferromagnetic Ising spin order at even larger values of h. Due to the large single-particle gaps for U > 6, we exclude residual (charge) fluctuations in the fermionic model to qualitatively modify this basic physics. We observe in Fig. 2 that the coexistence regime widens little with increasing U . This thinness of the coexistence regime within the considered parameter regime explains, why it was not resolved by the fermionic QMC methods [30]. It should however be noted, that within the effective quantum spin model we cannot identify the multicritical point of the underlying fermionic theory, which was located at (U, 1/h) ≈ (4.2, 0.28) in Ref. 30. Indeed, the multicritical IFP is of genuine fermionic nature, cf. the FRG results.
Discussion. We studied (2 + 1)-dimensional Dirac fermions coupled to two compatible OPs with O(N 1 ) and O(N 2 ) symmetry in the vicinity of the multicritical isotropic fixed point, providing an emergent O(N 1 + N 2 ) symmetry. Our FRG study predicts a strong irrelevance of perturbations of the O(N 1 + N 2 ) symmetry, which is consistent with the numerical result of an emergent symmetry in the related Dirac fermion lattice model [30]. Contrary to this numerical study, however, we do not find a single line of direct, continuous order-to-order transitions. Instead, we identify a robust region of coexistence of both orders, which is separated by continuous transitions from the other phases. In the case of the O(3) ⊕ Z 2 symmetry, relevant for Ref. 30, we could furthermore support this field-theoretical analysis by large-scale QMC simulations. Such a combined approach should be fruitful for uncovering the nature of quantum critical fermions coupled to bosonic fields under a wide range of conditions also in related systems.
where Γ (2) k is the Hessian in field space and R k is an infrared cutoff function, described in detail below. We employ the notation t = ln(k/Λ) for the renormalization group time. The graded trace STr involves a sum over continuous degrees of freedom as well as spinor indices, the latter entering with a minus sign. This exact evolution equation can be approximated by employing an ansatz based on the original form of the microscopic action, i.e., where we use the summation convention on repeated indices. This truncation is known as the extended local potential approximation or LPA . Furthermore we assume that the bosonic potential V k (φ, χ) is an analytic function of the O(N 1 ) and O(N 2 ) invariant quantities Anticipating the splitting of the bosonic degrees of freedom into longitudinal and transversal, we denote the combined independent bosonic fields as Φ := {φ L , φ T , χ L , χ T }. Moreover, we make use of the optimized linear regulator functions for both, bosons and fermions [2][3][4], i.e., where and Θ(x) is a step function. With these choices, we then feed the Ansatz of Eq. (2) into the Wetterich equation Eq. (1). The respective projections are described below, where to alleviate the notation, and unlike in the main text, we denote dimensionful quantities with a tilde.

(8)
The flow equation for the potential takes the form with threshold functions The potential u can be expanded in different ways corresponding to three different scenarios depending on whether the order-parameter fields are in the symmetric (SYM) or in the spontaneously symmetry broken regime (SSB): (1) SYM-SYM: Both OP fields remain in their symmetric phases, i.e., the minimum of the potential is at κ φ/χ = 0. In this case, the potential can be expanded as (2) SSB-SYM: The OP φ acquires an expectation value, while the other OP χ remains in its symmetric phase.
The related situation with the reversed role of the two OPs is referred to as SYM-SSB.

B. Projections on couplings
Equations for the bosonic couplings are obtained from Eq. (9) by taking the appropriate derivatives. In the SYM-SYM phase, we have In the SSB-SYM case, we have and finally, in the SSB-SSB case, we need to calculate In the latter case, the location of the minimum of the potential flows according to

C. Yukawa couplings
To be able to compare with results from the perturbative -expansion results, we take a projection along the longitudinal directions of the OPs. The coordinates are chosen such that φ 1 and χ 1 are the massive directions, this means The full flow equations thus obtained are with threshold functions given by where the index i, j ∈ {1, 2} referring to the two distinct Yukawa couplings is summed over when repeated and α ∈ {L, T } refers to the longitudinal or transverse components, respectively.

D. Anomalous dimensions
Finally, the set of equations is closed by considering the anomalous dimensions. A projection along the longitudinal components leads to , from which one gets with and the situation is qualitatively the same for all other combinations of N 1 and N 2 with N = N 1 + N 2 ≤ 5. Negative values of m 2 i , i ∈ {1, 2} correspond to a phase where the O(N i ) symmetry is spontaneously broken and thus the flow diagram Fig. 1 suggests a tendency towards phases where a single one of the symmetries spontaneously broken and the other one remains unbroken. This is represented by flows where m 2 1 < 0 stay in the m 2 2 > 0 region or vice versa, as indicated by the blue (green) stream lines. This preliminary analysis is incomplete for several reasons. First, the flows for m 2 i,k < 0 need to take into account that the minimum of the bosonic potential lies away from the origin in field space, a fact which is neglected in Eq. (11). Additionally, these streams can be misleading in the sense that there appears to be a saddle point for m 2 1 = m 2 2 < 0. Upon closer examination this is, in fact, not a fixed point of the whole set of flow equations.
To expand on this preliminary analysis, one needs to integrate the full FRG flows. We show a paradigmatic FRG flow in Fig. 2. In the upper panel, where the flow evolves through several regimes, as described above, and we switch the parametrization of the bosonic potential accordingly. We note that there is still a residual RG running in the deep infrared due to the presence of Goldstone modes in the SSB regimes, see also, for example, Refs. 5-7. For comparison, we also show the flow of the quartic couplings in the lower panel of Fig. 2.
Identifying the low-energy phase of the system for a particular choice of the initial conditions requires determining the regime of the effective potential which is adopted in the deep infrared, i.e., for k → 0, or equivalently, t → −∞.
To that end, we follow the values of the expectation values of the OPs, e.g., lim k→0 κ φ and lim k→0 κ χ .
In Fig. 2 both expectation values κ i are finite in the deep infrared and the flow ends up in the coexistence phase, where both O(N i ) symmetries are spontaneously broken. Similar flow trajectories like the one presented in Fig. 2 are found for other initial conditions and by systematically scanning parameter space in the vicinity of the IFP, we confirm that a system described by the effective theory of Eq. (2) indeed realizes all three possibilities mentioned around Eqs. (11)- (13). The phase where lim k→0 κ φ,k = lim k→0 κ χ,k = 0, i.e., where both symmetries are left unbroken, is connected to the Dirac semimetallic (DSM) phase. This phase can be further characterized in LPA∞ , i.e., to any order in the LPA , by nonvanishing universal fixed point values for all running couplings [6]. This means, in particular, that the Yukawa and quartic couplings satisfy g 2 1 = g 2 2 =: g 2 /v D and λ 2,0 = λ 1,1 = λ 0,2 =: λ/v D , and they flow to the universal values We note here that this just implies that the dimensionful quantities flow according to their canonical scaling, and the physically relevant quantity controlling the strength of interactions flows to zero in the infrared [8], Further, for a range of initial conditions, only one of the OPs acquires a finite expectation value. In the following discussion we fix the OP with the broken symmetry to be that corresponding to N 1 , i.e., φ, so that our statement can be rephrased as lim k→0 κ φ,k = 0 and lim k→0 κ χ,k = 0 .
This OP can thus be described as flowing to a Nambu-Goldstone infrared fixed point (of that symmetry) where the dimensionless expectation value satisfies while the other OP decouples and flows to a semimetallic-like fixed point (of the other symmetry), i.e, one where Eq. (30) holds. Finally in the vicinity of the line of exact O(N ) symmetry, both OPs acquire a finite expectation value. In particular, this means that crossing from an O(3)-broken phase into a Z 2 -broken phase takes place as a sequence of transitions in which both symmetries are broken after the first continuous transition. This is consistent with the observation in Ref. [9] that the single particle gap at the Dirac points, ∆ sp := 2(g 2 1 ρ φ + g 2 2 ρ χ ) does not close during the process and changes smoothly across such a transition. In the SSB regimes, we plot m 2 1 := ∆ φ and m 2 2 := ∆χ. Flows are initialized at t = 0 in the symmetric regime close to the IFP. At around t ≈ −9, the boson mass term m 2 1 → 0 and we switch to the symmetry-broken parametrization where κ φ = 0. Then at t ≈ −11.5 also m 2 2 → 0 and we switch to the parametrization of the boson potential where κ φ and κχ = 0. Lower panel: RG evolution for the dimensionless quartic couplings exhibiting NG plateaus in λ2,0 and λ1,1.
The phase associated to the coexistence region corresponds to an infrared behavior of the system that differs from a situation where all couplings flow together to an infrared fixed point of broken O(N 1 + N 2 ) symmetry. Moreover, the two OPs do not decouple inside the coexistence fan (as they do in the corresponding regions with only one of the symmetries broken) and this leads to an approximate O(N 1 + N 2 ) symmetry. Correspondingly, a snapshot of the joint (normalized) probability density P (φ a , χ b ) for the OPs across the coexistence region does not produce an exact circular pattern but an ellipse. Within our analysis, this can be confirmed from the behavior of the quantity where V k is the running bosonic potential, and whose level sets coincide with those of P (φ a , χ b ). A conour plot of f is shown in Fig. 3. The vicinity of the maximum at φ, χ = 0 has rotational symmetry between the order parameters and gives way to an elliptical shape away from it.

A. Numerical instability
The presence of Goldstone modes in phases with broken O(N i )-symmetry (N i > 1) leads to the limit lim k→0 κ i,k being ill-defined numerically close to the line of exact enlarged symmetry. Close to this line, the Goldstone modes force the expectation values of the OPs to switch constantly between zero and nonzero values along the flow in a manner that depends on the finite value of k at which the flow is stopped. That this is indeed the cause of the instability can be confirmed by "turning off" the Goldstone modes, i.e., by considering a model of two coupled Ising order parameters. In this case the expectation values of both order parameters converge to some nonzero values independent of the scale at which the flow is stopped, meeting continuously at the line of exact O(2) symmetry, as seen in Fig. 4. The residual nonmonotonic behavior as a function of couplings is nonuniversal and depends on the particular path chosen in parameter space.

III. DERIVATION OF THE EFFECTIVE QUANTUM SPIN MODEL
The Hamiltonian of the original fermion model from Ref. [9] reads